
C    COMBINED ROTOR-BOUNDARY LAYER PROGRAM

      COMMON/CODE/KODE,DELSRL,DELSRU
      COMMON/CTOBL/ALPH1,SPA,XMAX,TE,NTURBU,NTURBL,CTHETU,CTHETL
      COMMON/LNK1/XOML(100),YOML(100),SML(100),ETAML(100),NL,XOMU(100),
     1 YOMU(100),SMU(100),ETAMU(100),NU
      COMMON/LNK2/SMIN,SMOUT,BETAT,CSTAR,SIGMAO,GSTARO,GSTARI,GAMMA,
     * CONVER,RECONV
      COMMON/ROT1/V(4),BETAN,NKJ(5),RR(4)
      COMMON/C1/DUM1(7),NTURB,DUM2(6),CTHET,DUM3(910)

      DIMENSION X(100),Y(100)

    1 BETAP = 0.
      BETAM = 0.
C    CALCULATE UNROTATED BLADE COORDINATES
      CALL ROTORU
      CALL INPUT1
      CALB = 0.
C    CALCULATE ROTATED BLADE COORDINATES
    2 CALL ROTORR (CALB)
      CALB = 1.
C    NORMALIZE ROTATED BLADE COORDINATES TO CHORD LENGTH
      DO 3 I=1,NU
      X(I) = XOMU(I)/CSTAR
    3 Y(I) = YOMU(I)/CSTAR
      KODE = 0
      NTURB = NTURBU
      CTHET = CTHETU
      CALL INPUT2 (X,Y,SMU,ETAMU,NU)
C    CALCULATE UPPER SURFACE BOUNDARY LAYER
      CALL BLAYR

C    NORMALIZE ROTATED BLADE COORDINATES TO CHORD LENGTH
      DO 4 I=1,NL
      X(I) = XOML(I)/CSTAR
    4 Y(I) = YOML(I)/CSTAR
      KODE = 1
      NTURB = NTURBL
      CTHET = CTHETL
      CALL INPUT2 (X,Y,SML,ETAML,NL)
C    CALCULATE LOWER SURFACE BOUNDARY LAYER
      CALL BLAYR

C    CHECK FOR EQUAL SPACING AT INLET AND OUTLET
      DELTOT = (DELSRL + DELSRU)*CSTAR
      OUTSPA = GSTARO + DELTOT/COS(BETAT*CONVER)
      ELIMIT = .0001*CSTAR/XMAX
      IF (ABS(GSTARI - OUTSPA) .LE. ELIMIT) GO TO 5
      IF (OUTSPA .GE. GSTARI) BETAP = BETAT
      IF (OUTSPA .LT. GSTARI) BETAM = BETAT
      GSTARO = OUTSPA
      IF (BETAP .NE. 0.  .AND.  BETAM .NE. 0.) GO TO 6
C    CALCULATE NEW OUTLET FLOW ANGLE
      BETAT = -ARCOS(COS(BETAN)*SMIN/SMOUT*((1. + (GAMMA-1.)/2.*
     * SMOUT*SMOUT)/(1. + (GAMMA-1.)/2.*SMIN*SMIN))**((GAMMA+1.)/(2.*(
     * GAMMA-1.))) + DELTOT/GSTARI)*RECONV
      GO TO 2
    6 BETAT = (BETAP + BETAM)/2.
      GO TO 2

    5 WRITE (6,100) BETAT,GSTARI,OUTSPA
  100 FORMAT (//5X,19HITERATION COMPLETED/7X,11HBETA(OUT) =,F10.4,4H DEG
     *,10X,8HG*(IN) =,F10.5,10X,22HG*(OUT) + DEL(TOTAL) =,F10.5)
      CALB = 2.
      CALL ROTORR (CALB)
      GO TO 1

      END
$IBFTC INPU1
      SUBROUTINE INPUT1

C    BOUNDARY LAYER INPUT

      COMMON/LTU/LDAT
      COMMON/LNK2/UPMAC,SMOUT,ALP1M,CSTAR,SIGMAO,GSTARO,GSTARI,GAMMA,
     * CONVER,RECONV
      COMMON/CTOBL/ALPH1,SPA,XMAX,TE,NTURBU,NTURBL,CTHETU,CTHETL
      COMMON/C1/GAM,R,PTZ,TTZ,UPMACH,NST,NVP,NTURB,KPVM,KEM,KSMTH,
     1KSPLN,KLE,KATCH,CTHET,DLAM,TLAM,DTURB,TTURB,KPRE,KGRAD,KSDE,KLAM,
     2KMAIN,KPROF,X(100),Y(100),PRES(100),UE(100),ME(100),POPTZ(100),
     3VOVCR(100),TWAL(100),ETA(100)

      READ (LDAT,100) NTURBU,NTURBL,NVP,R,PTZ,TTZ,XMAX,TE,CTHETU,CTHETL
  100 FORMAT (3I3,1X,7(F9.5,1X))
      READ (LDAT,1010) KPRE,KGRAD,KSDE,KLAM,KMAIN,KPROF,KEM
 1010 FORMAT (7I5)

      GAM = GAMMA
      UPMACH = UPMAC

      RETURN
      END
$IBFTC ROTU
      SUBROUTINE ROTORU

C****CALCULATE THE UNROTATED BLADE COORDINATES

      COMMON/LTU/LDAT
      COMMON/FACTOR/PERM,GAMM1,GAMP1,GAM
      COMMON/ROT1/VIN,VOUT,VLOW,VUP,BETAN,NPER,KMAXN,KMAXO,JMAXN,JMAXO,
     1 RIN,ROUT,RLOW,RUP
      COMMON/ROT2/XLOW(400),YLOW(400),STML(400),ETALOW(400),KNDEX,
     1 XUP(400),YUP(400),STMU(400),ETAUP(400),JNDEX
      COMMON/LNK2/SMIN,SMOUT,BETAT,CSTAR,SIGMAO,GSTARO,GSTARI,GAMMA,
     * CONVER,RECONV

      LOGICAL ANGLE

      EXTERNAL FOFRS

      DATA LDAT/5/

      F(V,FN) = 2.*V - HALFPI*(PERM-1.) - 2.*(FN-1.)*DELV

    1 READ (LDAT,11) BETAN,VIN,VLOW,VUP,VOUT,BETAT,DELV,GAM
   11 FORMAT (8(F6.2,2X))

CC   CONVERSION FACTORS  AND  CONSTANTS
      GAMMA = GAM
      CONVER = .174532925E-01
      RECONV = 57.2957796
      HALFPI = 3.14159265/2.
C         ONE POINT WILL BE PRINTED FOR EVERY NPER POINTS CALCULATED
      NPER = 10
      IF (DELV .GE. 0.2) NPER = 1
      GAMP1 = (GAM + 1.)/2.
      GAMM1 = (GAM-1.)/2.
      PERM = SQRT(GAMP1/GAMM1)
      X0 = 1./PERM
      X2 = 0.999999999
      XINTL = (X0 + X2)/2.

      ANGLE = .TRUE.
      IF (VLOW .LE. AMIN1(VUP,VIN,VOUT)) GO TO 120
      WRITE (6,119)
  119 FORMAT (//31X,70HV(LOW) MUST BE LESS THAN OR EQUAL TO THE MINIMUM
     1OF V(UP),V(IN),V(OUT))
      ANGLE = .FALSE.

  120 IF (VUP .GE. AMAX1(VIN,VOUT)) GO TO 118
      WRITE (6,117)
  117 FORMAT (//33X,66HV(UP) MUST BE GREATER THAN OR EQUAL TO THE MAXIMU
     1M OF V(IN),V(OUT))
      ANGLE = .FALSE.

  118 VUMAX = HALFPI*(PERM-1.)*RECONV
      IF (VUP .LE. VUMAX) GO TO 116
      WRITE (6,115) VUMAX
  115 FORMAT (//41X,37HV(UP) MUST BE LESS THAN V(UP)(MAX) = ,F9.4,4H DEG
     1)
      ANGLE = .FALSE.

  116 IF (.NOT. ANGLE) GO TO 1

      WRITE (6,97)
   97 FORMAT (1H1,57X,17HDESIGN PARAMETERS)

CC   MISCELLANEOUS  CALCULATIONS
      DELV = DELV*CONVER
      FN = 1.
      V = VIN*CONVER
      DO 4 I=1,2
      FOFX = F(V,FN)
      CALL ROOT (X0,X2,XINTL,FOFX,FOFRS,X1)
      IF (I .EQ. 2) GO TO 4
      RIN = X1
      V = VOUT*CONVER
    4 CONTINUE
      ROUT = X1

      SSMIN = 1./RIN
      SSMOUT = 1./ROUT

      SMS = SSMIN
      I = 1
   16 SM = SQRT(((1./GAMP1)*SMS*SMS)/(1.-(GAMM1/GAMP1)*SMS*SMS))
      GO TO (17,18,19,20),I
   17 SMIN = SM
      SMS = SSMOUT
      I = 2
      GO TO 16
   18 SMOUT = SM
      DELV = DELV*RECONV

CC   PRINT  ALL  DESIGN  PARAMETERS
      WRITE (6,95) BETAN,VIN,VUP,VOUT,BETAT
   95 FORMAT (/2X,10HBETA(IN) =,F8.4,4H DEG,4X,7HV(IN) =,F8.4,4H DEG,6X,
     1 7HV(UP) =,F9.4,4H DEG,7X,8HV(OUT) =,F8.4,4H DEG,4X,11HBETA(OUT) =
     2,F9.4,4H DEG)
      WRITE (6,94) DELV, VLOW, GAM
   94 FORMAT (/20X,9HDELTA V =,F8.4,4H DEG,11X,8HV(LOW) =,F8.4,4H DEG,
     1 11X,7HGAMMA =,F8.4)

CC   CONVERT  FROM  DEGREES  TO  RADIANS
      VIN = VIN*CONVER
      VOUT = VOUT*CONVER
      VUP = VUP*CONVER
      VLOW = VLOW*CONVER
      BETAN = BETAN*CONVER
      DELV = DELV*CONVER

CC   CHOOSE LONGEST TRANSITION ARC OF LOWER SURFACE
      VNL = VIN - VLOW
      KMAXN = (VNL/DELV) + 0.5
      VOL = VOUT - VLOW
      KMAXO = (VOL/DELV) + 0.5
      KMN = MAX0(KMAXN,KMAXO)
      V = AMAX1(VIN,VOUT)

CC   CALCULATE  R*(LOW)=RLOW,  M*(LOW)=SSMLOW,  M(LOW)=SMLOW
      IF (VLOW .EQ. 0.0) GO TO 2
      FN = KMN + 1
      FOFX = F(V,FN)
      CALL ROOT (X0,X2,XINTL,FOFX,FOFRS,RLOW)
      GO TO 3
    2 RLOW = 1.0
    3 SSMLOW = 1./RLOW
      SMS = SSMLOW
      I = 3
      GO TO 16
   19 SMLOW = SM

CC   SET INITIAL POINTS FOR LOWER ARC CALCULATIONS
      KNDEX = KMN/NPER
      KDEX = KNDEX
      STML(KDEX+1) = SSMLOW
      XLOW(KDEX+1) = 0.0
      YLOW(KDEX+1) = RLOW
      ETALOW(KDEX+1) = 0.
      PHIKP1 = -(V-VLOW) + FLOAT(KMN)*DELV
      UMKP1 = ARSIN(SQRT(GAMP1*RLOW*RLOW - GAMM1))
      TXLO = XLOW(KDEX+1)
      TYLO = YLOW(KDEX+1)

CC   CHOOSE LONGEST TRANSITION ARC OF UPPER SURFACE
      VUT = VUP - VOUT
      JMAXO = (VUT/DELV)+0.5
      VUI = VUP - VIN
      JMAXN = (VUI/DELV)+0.5
      JMN = MAX0(JMAXO,JMAXN)
      V = AMIN1(VOUT,VIN)

CC   CALCULATE  R*(UP)=RUP,  M*(UP)=SSMUP,  M(UP)=MUP
      FN = -(JMN+1) + 2
      FOFX = F(V,FN)
      CALL ROOT (X0,X2,XINTL,FOFX,FOFRS,RUP)
      SSMUP = 1./RUP
      SMS = SSMUP
      I = 4
      GO TO 16
   20 SMUP = SM

CC   SET INITIAL POINTS FOR UPPER ARC CALCULATIONS
      JNDEX = JMN/NPER
      JDEX = JNDEX
      STMU(JDEX+1) = SSMUP
      XUP(JDEX+1) = 0.0
      YUP(JDEX+1) = RUP
      ETAUP(JDEX+1) = 0.
      PHIJP1 = -(VUP-V) + FLOAT(JMN)*DELV
      UMJP1 = ARSIN(SQRT(GAMP1*RUP*RUP - GAMM1))
      TXUP = XUP(JDEX+1)
      TYUP = YUP(JDEX+1)

      IF (VIN .EQ. VLOW  .AND.  VLOW .EQ. VOUT) GO TO 100

C****CALCULATE COORDINATES FOR LOWER TRANSITION ARC - UNROTATED
      KDEX = KNDEX + 1
      NUM = 0
      V = AMAX1(VIN,VOUT)
      KEEP = KMN + 1
      DO 30 KK=1,KMN
      K = KEEP - KK
      NUM = NUM + 1
      PHIK = PHIKP1 - DELV
      FN = K
      FOFX = F(V,FN)
      CALL ROOT (X0,X2,XINTL,FOFX,FOFRS,TR)
      TX = TR*SIN(PHIK)
      TY = TR*COS(PHIK)
      EMWK = TAN(-PHIKP1)
      UMK = ARSIN(SQRT(GAMP1*TR*TR - GAMM1))
      EMK = -TAN((PHIK+UMK+PHIKP1+UMKP1)/2.)
      TEMP = TYLO - EMWK*TXLO
      TEMPP = TY - EMK*TX
      TEMPPP = EMK - EMWK
      TXLO = (TEMP - TEMPP)/TEMPPP
      TYLO = ((EMK*TEMP) - (EMWK*TEMPP))/TEMPPP
      PHIKP1 = PHIK
      UMKP1 = UMK

CC        SAVE EVERY "NPER-TH" POINT
      N = NUM - (NUM/NPER)*NPER
      IF (N .GT. 0) GO TO 30
      KDEX = KDEX - 1
      STML(KDEX) = 1./TR
      XLOW(KDEX) = TXLO
      YLOW(KDEX) = TYLO
      ETALOW(KDEX) = -PHIK

   30 CONTINUE

  100 IF (VIN .EQ. VUP  .AND.  VUP .EQ. VOUT) GO TO 200

C****CALCULATE COORDINATES FOR UPPER TRANSITION ARC - UNROTATED
      JDEX = JNDEX + 1
      NUM = 0
      V = AMIN1(VOUT,VIN)
      JEEP = JMN + 1
      DO 41 JJ=1,JMN
      J = JEEP - JJ
      NUM = NUM + 1
      PHIJ = PHIJP1 - DELV
      FN = -J + 2
      FOFX = F(V,FN)
      CALL ROOT (X0,X2,XINTL,FOFX,FOFRS,TR)
      TX = TR*SIN(PHIJ)
      TY = TR*COS(PHIJ)
      EMWJ = TAN(-PHIJP1)
      UMJ = ARSIN(SQRT(GAMP1*TR*TR - GAMM1))
      EMJ = TAN((-PHIJ+UMJ-PHIJP1+UMJP1)/2.)
      TEMP = TYUP - EMWJ*TXUP
      TEMPP = TY - EMJ*TX
      TEMPPP = EMJ - EMWJ
      TXUP = (TEMP - TEMPP)/TEMPPP
      TYUP = ((EMJ*TEMP) - (EMWJ*TEMPP))/TEMPPP
      PHIJP1 = PHIJ
      UMJP1 = UMJ

CC        SAVE EVERY "NPER-TH" POINT
      N = NUM - (NUM/NPER)*NPER
      IF (N .GT. 0) GO TO 41
      JDEX = JDEX - 1
      STMU(JDEX) = 1./TR
      XUP(JDEX) = TXUP
      YUP(JDEX) = TYUP
      ETAUP(JDEX) = -PHIJ

   41 CONTINUE

CC   MISCELLANEOUS OUTPUT
  200 WRITE (6,622)
  622 FORMAT (//54X,24HMISCELLANEOUS PARAMETERS//)
      WRITE (6,1000) SSMIN,SMIN,SMOUT,SSMOUT
 1000 FORMAT (/25X,8HM*(IN) =,F9.4,3X,7HM(IN) =,F9.4,10X,8HM(OUT) =,F9.4
     1,5X,9HM*(OUT) =,F9.4)
      WRITE (6,1001) SSMLOW,SMLOW,SMUP,SSMUP
 1001 FORMAT (/25X,9HM*(LOW) =,F9.4,2X,8HM(LOW) =,F9.4,11X,7HM(UP) =,
     * F9.4,6X,8HM*(UP) =,F9.4)

      KNDEX = KNDEX + 1
      JNDEX = JNDEX + 1

      RETURN
      END
$IBFTC ROO
      SUBROUTINE ROOT (X0,X2,XINTL,FOFX,FUNC,X1)

      DOUBLE PRECISION X,XX0,XX2

C    WE ARE SEEKING AN X SUCH THAT FUNC(X) = FOFX WHERE FOFX IS A KNOWN
C    FUNCTIONAL VALUE
C      1  LOCATE FOFX IN (F0,FX) OR (FX,F2) WHERE FX IS THE PREVIOUS
C         APPROXIMATION TO FOFX
C      2  LET X = 1/2(XX0+X) OR X = 1/2(X+XX2)
C      3  IS FUNC(X) = FOFX ?  IF NOT, REPEAT PROCEDURE

      XX0 = X0
      XX2 = X2
      F0 = FUNC(XX0)
      F2 = FUNC(XX2)
      IF ( FOFX .LT. F0 .AND. FOFX .LT. F2  .OR.  FOFX .GT. F0 .AND.
     1FOFX .GT. F2 ) GO TO 1005

      IF (ABS(FOFX-F0) .LE. .00001) GO TO 1007
      IF (ABS(FOFX-F2) .LE. .00001) GO TO 1008

      X = XINTL
      KOUNT = 0
 1000 X1 = X
      KOUNT = KOUNT + 1
      A = FOFX - F2
      FX = FUNC(X)
      IF (KOUNT .GE. 60) WRITE (6,1004) KOUNT,X,FX,FOFX
 1004 FORMAT (1HL,9H  KOUNT  ,G16.9,9H    X    ,G16.9,9H   FX    ,G16.9,
     19H   FOFX  ,G16.9)
      IF (ABS(FX-FOFX) .LE. .00001) RETURN
      IF (KOUNT .EQ. 75) GO TO 1002
      IF (A*(FX-FOFX) .LT. 0.) GO TO 1001
      XX0 = X
      X = (X+XX2)/2.
      GO TO 1000
 1001 XX2 = X
      X = (XX0+X)/2.
      F2 = FX
      GO TO 1000

 1002 WRITE (6,1003)
 1003 FORMAT (//30X,62H75 ITERATIONS HAVE BEEN PERFORMED WITHOUT CONVERG
     1ING TO A ROOT)
      RETURN

 1005 WRITE (6,1006) FOFX
 1006 FORMAT (//10X,7HF(X) = ,G16.9,31H IS OUTSIDE OF SPECIFIED LIMITS)
      RETURN

 1007 X1 = X0
      RETURN

 1008 X1 = X2
      RETURN

      END
$IBFTC FELI
      FUNCTION FOFRS (X)

      DOUBLE PRECISION X
      COMMON/FACTOR/PERM,GAMM1,GAMP1,GAM

      ARG1 = 2.*GAMM1/(X*X) - GAM
      ARG2 = 2.*GAMP1*X*X - GAM
      IF (ABS(ARG1) .GT. 1.0  .OR.  ABS(ARG2) .GT. 1.0) WRITE (6,1) ARG1
     1,ARG2
    1 FORMAT (//14X,61HARGUMENT OF ARCSIN IS OUTSIDE DOMAIN OF DEFINITIO
     1N    ARG1 = ,G16.9,11H    ARG2 = ,G16.9)

      FOFRS = PERM*ARSIN(ARG1) + ARSIN(ARG2)

      RETURN
      END
$IBFTC ROTR
      SUBROUTINE ROTORR (CALB)

C    CALCULATE  (1) ROTATED BLADE COORDINATES
C               (2) INLET AND OUTLET BLADE SPACING

      DIMENSION XLOWN(80),YLOWN(80),SMLN(80),ETALN(80),XCLOW(35),
     1 YCLOW(35),ETACL(35),XLOWO(80),YLOWO(80),SMLO(80),ETALO(80)
      DIMENSION XSN(11),YSN(11),ETASN(11),XUPN(80),YUPN(80),SMUN(80),
     1 ETAUN(80),XCUP(35),YCUP(35),ETACU(35),XUPO(80),YUPO(80),SMUO(80),
     2 ETAUO(80),XSO(11),YSO(11),ETASO(11)

      COMMON/LNK1/XOML(100),YOML(100),SML(100),ETAML(100),NL,XOMU(100),
     1 YOMU(100),SMU(100),ETAMU(100),NU
      COMMON/ROT1/VIN,VOUT,VLOW,VUP,BETAN,NPER,KMAXN,KMAXO,JMAXN,JMAXO,
     1 RIN,ROUT,RLOW,RUP
      COMMON/ROT2/XLOW(400),YLOW(400),STML(400),ETALOW(400),KNDEX,
     1 XUP(400),YUP(400),STMU(400),ETAUP(400),JNDEX
      COMMON/LNK2/SMIN,SMOUT,BETAT,CSTAR,SIGMAO,GSTARO,GSTARI,GAMMA,
     * CONVER,RECONV

      IF (CALB .NE. 0.) GO TO 2
      WRITE (6,100) BETAT
  100 FORMAT (//10X,15HINPUT BETA(0) =,F10.4,4H DEG///)
      GO TO 3

    2 IF (CALB .NE. 1.) GO TO 3
      WRITE (6,101) BETAT,GSTARO
  101 FORMAT (//10X,20HCALCULATED BETA(0) =,F10.4,4H DEG,5X,14HWHEN G*(O
     1UT) =,F10.5)

    3 BETAT = BETAT*CONVER

      ALPHLN = (VIN - VLOW) - BETAN
      ALPHLO = -(VOUT - VLOW) - BETAT
      IF (ALPHLN .LE. 0.  .AND.  ALPHLO .GE. 0.) GO TO 4
      WRITE (6,102)
  102 FORMAT (//27X,79HV(LOW) MUST BE GREATER THAN OR EQUAL TO V(IN) - B
     1ETA(IN) AND V(OUT) + BETA(OUT))
      STOP

C****CALCULATE COORDINATES FOR LOWER TRANSITION ARC - ROTATED
    4 SINALN = SIN(ALPHLN)
      COSALN = COS(ALPHLN)
      ABSALN = ABS(ALPHLN)
      SINALO = SIN(ALPHLO)
      COSALO = COS(ALPHLO)
      ABSALO = ABS(ALPHLO)
      KDEX = KNDEX + 1
      KN = (KMAXN/NPER) + 2
      KO = (KMAXO/NPER) + 2

      DO 5 KK=1,KNDEX
      K = KDEX - KK
      KN = KN - 1
      KO = KO - 1
      IF (KN .LE. 0) GO TO 6
      SMLN(KN) = STML(K)
      XLOWN(KN) = YLOW(K)*SINALN + XLOW(K)*COSALN
      YLOWN(KN) = YLOW(K)*COSALN - XLOW(K)*SINALN
      ETALN(KN) = ETALOW(K) + ABSALN
    6 IF (KO .LE. 0) GO TO 5
      SMLO(KO) = STML(K)
      XLOWO(KO) = YLOW(K)*SINALO - XLOW(K)*COSALO
      YLOWO(KO) = YLOW(K)*COSALO + XLOW(K)*SINALO
      ETALO(KO) = ETALOW(K) + ABSALO
    5 CONTINUE

      ALPHUN = (VUP - VIN) - BETAN
      ALPHUO = -(VUP - VOUT) - BETAT
      IF (ALPHUN .LE. 0.  .AND.  ALPHUO .GE. 0.) GO TO 7
      WRITE (6,103)
  103 FORMAT (//2X,75HV(UP) MUST BE LESS THAN OR EQUAL TO V(IN) + BETA(I
     1N) AND V(OUT) - BETA(OUT))
      STOP

C****CALCULATE COORDINATES FOR UPPER TRANSITION ARC - ROTATED
    7 SINAUN = SIN(ALPHUN)
      COSAUN = COS(ALPHUN)
      ABSAUN = ABS(ALPHUN)
      SINAUO = SIN(ALPHUO)
      COSAUO = COS(ALPHUO)
      ABSAUO = ABS(ALPHUO)
      JDEX = JNDEX + 1
      JN = (JMAXN/NPER) + 2
      JO = (JMAXO/NPER) + 2

      DO 8 JJ=1,JNDEX
      J = JDEX - JJ
      JO = JO - 1
      JN = JN - 1
      IF (JO .LE. 0) GO TO 9
      SMUO(JO) = STMU(J)
      XUPO(JO) = YUP(J)*SINAUO - XUP(J)*COSAUO
      YUPO(JO) = YUP(J)*COSAUO + XUP(J)*SINAUO
      ETAUO(JO) = ETAUP(J) + ABSAUO
    9 IF (JN .LE. 0) GO TO 8
      SMUN(JN) = STMU(J)
      XUPN(JN) = YUP(J)*SINAUN + XUP(J)*COSAUN
      YUPN(JN) = YUP(J)*COSAUN - XUP(J)*SINAUN
      ETAUN(JN) = ETAUP(J) + ABSAUN
    8 CONTINUE

C****CALCULATE THE INLET AND OUTLET (DIMENSIONLESS) BLADE SPACING
      YLASTI = YUPN(1) + TAN(BETAN)*(XLOWN(1) - XUPN(1))
      GSTARI = YLOWN(1) - YLASTI
      YLASTO = YUPO(1) + TAN(BETAT)*(XLOWO(1) - XUPO(1))
      GSTARO = YLOWO(1) - YLASTO

      DALPH = 5.*CONVER
C****CIRCULAR ARC (LOWER)
      ALPH = ALPHLO + DALPH
      ALPLOW = ALPHLN
      KOUNT = 0
   10 KOUNT = KOUNT + 1
      XCLOW(KOUNT) = RLOW*SIN(ALPLOW)
      YCLOW(KOUNT) = RLOW*COS(ALPLOW)
      ETACL(KOUNT) = ABS(ALPLOW)
      ALPLOW = ALPLOW + DALPH
      IF (ABS(ALPH-ALPLOW) .LE. .001) GO TO 11
      IF (ALPHLO .LT. ALPLOW  .AND.  ALPLOW .LT. ALPH) ALPLOW = ALPHLO
      GO TO 10
   11 NCL = KOUNT

C****CIRCULAR ARC (UPPER)
      ALPH = ALPHUO + DALPH
      ALPHUP = ALPHUN
      KOUNT = 0
   12 KOUNT = KOUNT + 1
      XCUP(KOUNT) = RUP*SIN(ALPHUP)
      YCUP(KOUNT) = RUP*COS(ALPHUP)
      ETACU(KOUNT) = ABS(ALPHUP)
      ALPHUP = ALPHUP + DALPH
      IF (ABS(ALPH-ALPHUP) .LE. .001) GO TO 13
      IF (ALPHUO .LT. ALPHUP  .AND.  ALPHUP .LT. ALPH) ALPHUP = ALPHUO
      GO TO 12
   13 NCU = KOUNT

C****CALCULATE COORDINATES FOR STRAIGHT LINE PORTION OF UPPER ARC
      KOUNT = 1
      IF (XLOWN(1) .LE. XUPN(1)) GO TO 15
      WRITE (6,106)
  106 FORMAT (8H0(ROTRR),4X,71HUPPER SURFACE INLET LONGER THAN LOWER SUR
     *FACE INLET  -  CASE TERMINATED)
      STOP
   15 DELXI = (XUPN(1) - XLOWN(1))/10.
      IF (XLOWO(1) .GE. XUPO(1)) GO TO 16
      WRITE (6,105)
  105 FORMAT (8H0(ROTRR),4X,73HUPPER SURFACE OUTLET LONGER THAN LOWER SU
     *RFACE OUTLET  -  CASE TERMINATED)
      STOP
   16 DELXO = (XLOWO(1) - XUPO(1))/10.
      XSIN = XUPN(1)
      YSIN = YUPN(1)
      XSOUT = XUPO(1)
      YSOUT = YUPO(1)
      XSN(KOUNT) = XSIN
      YSN(KOUNT) = YSIN
      ETASN(KOUNT) = ETAUN(1)
      XSO(KOUNT) = XSOUT
      YSO(KOUNT) = YSOUT
      ETASO(KOUNT) = ETAUO(1)
      TANBN = TAN(BETAN)
      ABSBN = ABS(BETAN)
      TANBO = TAN(BETAT)
      ABSBT = ABS(BETAT)
      DO 14 I=2,11
      XSIN = XSIN - DELXI
      XSOUT = XSOUT + DELXO
      YSIN = YUPN(1) + TANBN*(XSIN - XUPN(1))
      YSOUT = YUPO(1) + TANBO*(XSOUT - XUPO(1))
      ETASN(I) = ABSBN
      XSN(I) = XSIN
      YSN(I) = YSIN
      ETASO(I) = ABSBT
      XSO(I) = XSOUT
   14 YSO(I) = YSOUT

C****PREPARE DATA FOR BOUNDARY LAYER PROGRAM
      CSTAR = SQRT((XLOWO(1) - XLOWN(1))**2 + (YLOWO(1) - YLOWN(1))**2)
      SIGMAI = CSTAR/GSTARI
      SIGMAO = CSTAR/GSTARO
      WRITE (6,104) GSTARI,SIGMAI,CSTAR,GSTARO,SIGMAO
  104 FORMAT (//39X,8HG*(IN) =,F9.4,37X,11HSIGMA(IN) =,F9.4/69X,4HC* =,
     * F9.4/39X,9HG*(OUT) =,F9.4,36X,12HSIGMA(OUT) =,F9.4//)
      BETAT = BETAT*RECONV
      SSMLOW = 1./RLOW
      SSMUP = 1./RUP
      SSMIN = 1./RIN
      SSMOUT = 1./ROUT
      IF (CALB .EQ. 2.) RETURN

C    STORE ROTATED BLADE COORDINATES FOR LOWER SURFACE
      KN = (KMAXN/NPER) + 1
      DO 310 I=1,KN
      XOML(I) = XLOWN(I)
      YOML(I) = YLOWN(I)
      SML(I) = SMLN(I)
      ETAML(I) = ETALN(I)
  310 CONTINUE
      NL = KN
      NCL = NCL-1
      DO 311 I=2,NCL
      NL = NL + 1
      XOML(NL) = XCLOW(I)
      YOML(NL) = YCLOW(I)
      SML(NL) = SSMLOW
      ETAML(NL) = ETACL(I)
  311 CONTINUE
      KO = (KMAXO/NPER) + 1
      DO 312 I=1,KO
      J = KO + 1 - I
      NL = NL + 1
      XOML(NL) = XLOWO(J)
      YOML(NL) = YLOWO(J)
      SML(NL) = SMLO(J)
      ETAML(NL) = ETALO(J)
  312 CONTINUE

C    STORE ROTATED BLADE COORDINATES FOR UPPER SURFACE
      DO 313 I=1,10
      J = 12 - I
      XOMU(I) = XSN(J)
      YOMU(I) = YSN(J)
      SMU(I) = SSMIN
      ETAMU(I) = ETASN(J)
  313 CONTINUE
      NU = 10
      JN = (JMAXN/NPER) + 1
      DO 314 I=1,JN
      NU = NU + 1
      XOMU(NU) = XUPN(I)
      YOMU(NU) = YUPN(I)
      SMU(NU) = SMUN(I)
      ETAMU(NU) = ETAUN(I)
  314 CONTINUE
      NCU = NCU-1
      DO 315 I=2,NCU
      NU = NU + 1
      XOMU(NU) = XCUP(I)
      YOMU(NU) = YCUP(I)
      SMU(NU) = SSMUP
      ETAMU(NU) = ETACU(I)
  315 CONTINUE
      JO = (JMAXO/NPER) + 1
      DO 316 I=1,JO
      J = JO + 1 - I
      NU = NU + 1
      XOMU(NU) = XUPO(J)
      YOMU(NU) = YUPO(J)
      SMU(NU) = SMUO(J)
      ETAMU(NU) = ETAUO(J)
  316 CONTINUE
      DO 317 I=2,11
      NU = NU + 1
      XOMU(NU) = XSO(I)
      YOMU(NU) = YSO(I)
      SMU(NU) = SSMOUT
      ETAMU(NU) = ETASO(I)
  317 CONTINUE

      RETURN
      END
$IBFTC INPU2
      SUBROUTINE INPUT2 (XM,YM,VVCR,ETAM,NN)

      LOGICAL ERROR,TRANS,SEPRN

      COMMON/CODE/KODE,DELSRL,DELSRU
      COMMON/LNK2/UPMAC,SMOUT,ALP1M,CSTAR,SIGMAO,GSTARO,GSTARI,GAMMA,
     * CONVER,RECONV
      COMMON/CTOBL/ALPH1,SPA,XMAX,TE,NTURBU,NTURBL,CTHETU,CTHETL
      COMMON/C1/GAM,R,PTZ,TTZ,UPMACH,NST,NVP,NTURB,KPVM,KEM,KSMTH,
     1KSPLN,KLE,KATCH,CTHET,DLAM,TLAM,DTURB,TTURB,KPRE,KGRAD,KSDE,KLAM,
     2KMAIN,KPROF,X(100),Y(100),PRES(100),UE(100),ME(100),POPTZ(100),
     3VOVCR(100),TWAL(100),ETA(100)
      COMMON/C3/XOM(100),YOM(100),S(100),SOL(100),AE(100),TSE(100),
     1TAWL(100),TAWT(100),TBAR(100),RW(100),SW(100),SUTHL(100),
     2RHSW(100),RHSE(100),HEADW(100),HEADE(100),NUW(100),MUBAR(100),
     3AA(100),BB(100),FF(100),DUDS(100),DMDS(100),DMDL(100)
      COMMON/C9/ERROR,TRANS,SEPRN

      DIMENSION XM(NN),YM(NN),VVCR(NN),ETAM(NN)

C   THE V/VCR DISTRIBUTION IS RECEIVED FROM THE ROTOR PROGRAM
C    IT WILL BE USED TO CALCULATE THE PRESSURE DISTRIBUTION

      ERROR = .FALSE.
      TRANS = .FALSE.
      SEPRN = .FALSE.
C    BOUNDARY LAYER SETUP
      KATCH = 0
      KSPLN = 1
      KLE = 1
      TLAM = 0.
      DLAM = 0.
      TTURB = 0.
      DTURB = 0.

      NST = NN
      DO 2 I=1,NST
      XOM(I) = XM(I)
      YOM(I) = YM(I)
      X(I) = XOM(I)*XMAX
      Y(I) = YOM(I)*XMAX
      TWAL(I) = TTZ
      ETA(I) = ETAM(I)
    2 VOVCR(I) = VVCR(I)
      ALPH1 = -ALP1M
      SPA = XMAX/SIGMAO

      IF (KODE .NE. 0) GO TO 3
      WRITE (6,100)
  100 FORMAT (1H1,53X,23H***  UPPER SURFACE  ***)
      GO TO 4
    3 WRITE (6,101)
  101 FORMAT (1H1,53X,23H***  LOWER SURFACE  ***)
    4 WRITE (6,1020) NST,NVP,GAM,R,PTZ,TTZ,UPMACH,XMAX,NTURB,CTHET,TE
 1020 FORMAT (1H0/5X,22HBOUNDARY LAYER - INPUT///5X,3HNST,5X,3HNVP,7X,
     * 3HGAM,11X,1HR,12X,3HPTZ,11X,3HTTZ,9X,6HUPMACH,8X,4HXMAX,7X,
     * 5HNTURB,6X,6HCTHETA,7X,2HTE/4X,2(I3,5X),F7.3,3(5X,F9.2),5X,F8.4,
     * 5X,F8.5,6X,I3,7X,F7.3,5X,F7.5///)
      WRITE (6,1090) KPRE,KGRAD,KSDE,KLAM,KMAIN,KPROF,KEM
 1090 FORMAT(/6X,4HKPRE,7X,5HKGRAD,7X,4HKSDE,8X,4HKLAM,7X,5HKMAIN,7X,5HK
     1PROF,7X,5H KEM /7X,I2,9X,I2,10X,I2,10X,I2,9X,I2,10X,I2,10X,I2///)
      WRITE (6,1032) (XOM(I),YOM(I),VVCR(I),TWAL(I),I=1,NST)
 1032 FORMAT(8X,3HXOM,7X,3HYOM,9X,4HVVCR,9X,4HTWAL/(5X,F8.5,2X,F8.5,2X,
     1F11.5,3X,F9.2))

      IF(NST.GT.100.OR.NTURB.GT.NST.OR.KEM.LT.0.OR.KEM.GT.1.OR.KSPLN.LT.
     10.OR.KSPLN.GT.1.OR.KLE.LT.0.OR.KLE.GT.1.OR.KATCH.LT.0.OR.KATCH.GT.
     21) GO TO 70
      RETURN

   70 ERROR = .TRUE.
      WRITE(6,1170)
 1170 FORMAT(////10X,48HERROR IN INPUT DATA.  RECHECK INPUT INSTRUCTIONS
     1)
      RETURN

      END
$IBFTC BLAYR
      SUBROUTINE BLAYR

C    BOUNDARY LAYER MAIN PROGRAM

      LOGICAL ERROR,TRANS,SEPRN
      REAL ME

      COMMON/CODE/KODE,DELSRL,DELSRU
      COMMON/CTOBL/AL,SP,XMAX,TE,NTURBU,NTURBL,CHTETU,CTHETL
      COMMON/C1/GAM,R,PTZ,TTZ,UPMACH,NST,NVP,NTURB,KPVM,KEM,KSMTH,
     1KSPLN,KLE,KATCH,CTHET,DLAM,TLAM,DTURB,TTURB,KPRE,KGRAD,KSDE,KLAM,
     2KMAIN,KPROF,X(100),Y(100),PRES(100),UE(100),ME(100),POPTZ(100),
     3VOVCR(100),TWAL(100),ETA(100)
      COMMON/C2/PSZ,TSZ,UZ,ASZ,ATZ,RHSZ,RHTZ,MUSZ,MUTZ,NUSZ,NUTZ,CP,
     1PR,TC,ARCL
      COMMON/C3/XOM(100),YOM(100),S(100),SOL(100),AE(100),TSE(100),
     1TAWL(100),TAWT(100),TBAR(100),RW(100),SW(100),SUTHL(100),
     2RHSW(100),RHSE(100),HEADW(100),HEADE(100),NUW(100),MUBAR(100),
     3AA(100),BB(100),FF(100),DUDS(100),DMDS(100),DMDL(100)
      COMMON/C4/THET(100),DELSR(100),DELTA(100),FORM(100),
     1FORMI(100),FORMTR(100),RTH(100),RHTI(100),CF(100),
     2TAUW(100),NUSS(100),DTDY(100),HTRAN(100),CRN(100)
      COMMON/C9/ERROR,TRANS,SEPRN

      CALL PRECAL
      IF (ERROR) RETURN
      CALL LAMNAR
      IF (ERROR) RETURN
      IF (SEPRN) GO TO 20
      IF (.NOT. TRANS) GO TO 20
      CALL TURBLN
      IF (ERROR) RETURN
   20 CALL PROFIL
      IF (KODE.EQ.1) GO TO 9
      DELS = DELSR(NST)
      THETS = THET(NST)
      KODE = 1
      RETURN
    9 DELP = DELSR(NST)
      THETP = THET(NST)
      ALPH1 = AL
      SPA = SP
      CALL AFMIX (ALPH1,DELS,DELP,THETS,THETP,TE,SPA,ME(NST))
      KODE = 0
      RETURN

      END
$IBFTC PRECL
      SUBROUTINE PRECAL

      COMMON/GAMPM/GMP1,GMM1
      COMMON/CTOBL/DUM(2),XMAX,DUM2(5)
      COMMON/C1/GAM,R,PTZ,TTZ,UPMACH,NST,NVP,NTURB,KPVM,KEM,KSMTH,
     1KSPLN,KLE,KATCH,CTHET,DLAM,TLAM,DTURB,TTURB,KPRE,KGRAD,KSDE,KLAM,
     2KMAIN,KPROF,X(100),Y(100),PRES(100),UE(100),ME(100),POPTZ(100),
     3VOVCR(100),TWAL(100),ETA(100)
      COMMON/C2/PSZ,TSZ,UZ,ASZ,ATZ,RHSZ,RHTZ,MUSZ,MUTZ,NUSZ,NUTZ,CP,
     1PR,TC,ARCL
      COMMON/C3/XOM(100),YOM(100),S(100),SOL(100),AE(100),TSE(100),
     1TAWL(100),TAWT(100),TBAR(100),RW(100),SW(100),SUTHL(100),
     2RHSW(100),RHSE(100),HEADW(100),HEADE(100),NUW(100),MUBAR(100),
     3AA(100),BB(100),FF(100),DUDS(100),DMDS(100),DMDL(100)
      COMMON/C9/ERROR,TRANS,SEPRN
      DIMENSION SDER(100),CMU(20),CPR(20),CTC(20)
      REAL MUSZ,MUTZ,NUSZ,NUTZ,MUSLE,MUSLM,ME,NUW,MUBAR
      LOGICAL ERROR,TRANS,SEPRN
C
C  READ DATA FOR MU, PR, AND TC CURVE FITS
C
      DATA(CMU(I),I=1,5)/-.01945170,1.3019531,-.34511323,
     1.068277826,-.00566593/
      DATA(CPR(I),I=1,5)/.8557,-.234136,.1078624,
     1-.0236214,.00202863/
      DATA(CTC(I),I=1,5)/-.03839323,1.2697427,-.30911252,
     1.08743781,-.009674725/
C
      GMP1 = GAM + 1.
      GMM1 = GAM - 1.

C  INITIALIZE STATIC AND TOTAL PARAMETERS
C
      TSLE= 518.688
      TSLM= 288.160
      MUSLE= 3.711402E-7
      MUSLM= 1.777029E-5
      TCSLE= 3.202206E-3
      TCSLM= 2.561796E-2
      TSZ = TTZ/(1.+GMM1/2.*UPMACH**2)
      PSZ = PTZ*(TSZ/TTZ)**(GAM/GMM1)
      RHSZ= PSZ/R/TSZ
      RHTZ= PTZ/R/TTZ
      ASZ= SQRT(GAM*R*TSZ)
      ATZ= SQRT(GAM*R*TTZ)
      UZ= UPMACH*ASZ
      CP = R*GAM/GMM1
      IF (KEM.EQ.1) GO TO 10
      TCON= 198.60
      TR1= TSZ/TSLE
      TR2= TTZ/TSLE
      GO TO 20
   10 TCON= 110.33
      TR1= TSZ/TSLM
      TR2= TTZ/TSLM
   20 CALL CURVFT(CPR,PR,TR1,0,4,0)
      CALL CURVFT(CTC,TC,TR1,0,4,0)
      CALL CURVFT(CMU,MUSZ,TR1,0,4,0)
      CALL CURVFT(CMU,MUTZ,TR2,0,4,0)
      IF (KEM.EQ.1) GO TO 30
      TC= TC*TCSLE
      MUSZ= MUSZ*MUSLE
      MUTZ= MUTZ*MUSLE
      GO TO 40
   30 TC= TC*TCSLM
      MUSZ= MUSZ*MUSLM
      MUTZ= MUTZ*MUSLM
   40 NUSZ= MUSZ/RHSZ
      NUTZ= MUTZ/RHTZ
C
C  CALCULATE GEOMETRY RATIOS AND ARC LENGTHS
C
      S(1)= 0.
      DO 50 I=2,NST
   50 S(I)= S(I-1)+SQRT((X(I)-X(I-1))**2+(Y(I)-Y(I-1))**2)
      ARCL= S(NST)
      DO 60 I=1,NST
   60 SOL(I)= S(I)/ARCL
C
C  CALCULATE PRES,UE,ME,POPTZ,AND VOVCR AT EACH STATION
C
C    VELOCITY OVER CRITICAL VELOCITY GIVEN AS INPUT
  150 DO 160 I=1,NST
      POPTZ(I) = (1.-GMM1/GMP1*VOVCR(I)**2)**(GAM/GMM1)
      UE(I) = SQRT(2.*GAM/GMM1*PTZ/RHTZ*(1.-POPTZ(I)**(GMM1/GAM)
     *))
      TSE(I) = TTZ-UE(I)**2/(2.*CP)
      AE(I) = SQRT(GAM*R*TSE(I))
      ME(I) = UE(I)/AE(I)
  160 PRES(I) = POPTZ(I)*PTZ
C
C  PRINT INITIAL CALCULATED PARAMETERS
C
  170 WRITE(6,1000)
      WRITE(6,1010) PSZ,TSZ,UZ,ASZ,ATZ,RHSZ,RHTZ,MUSZ,MUTZ,NUSZ,NUTZ,CP,
     1PR,TC,ARCL
      WRITE(6,1020) (I,PRES(I),UE(I),ME(I),POPTZ(I),VOVCR(I),I=1,NST)
C
C
C  PRINT GEOMETRY PARAMETERS
C
  200 IF (KPRE.NE.1) GO TO 210
      WRITE(6,1030) (I,X(I),Y(I),S(I),XOM(I),YOM(I),SOL(I),I=1,NST)
C
C  CALCULATE OTHER NECESSARY PARAMETERS AT EACH STATION
C
  210 DO 220 I=1,NST
      TEM1 = 1.+GMM1/2.*ME(I)**2
      RHSW(I)= PRES(I)/R/TWAL(I)
      RHSE(I)= PRES(I)/R/TSE(I)
      HEADW(I)= .5*RHSW(I)*UE(I)**2
      HEADE(I)= .5*RHSE(I)*UE(I)**2
      SW(I)= TWAL(I)/TTZ-1.
      SUTHL(I)= SQRT(TWAL(I)/TTZ)*(TTZ+TCON)/(TWAL(I)+TCON)
      NUW(I) = SUTHL(I)*NUTZ*(1.+SW(I))**2*TEM1**(GAM/GMM1)
      RW(I)= UE(I)*S(I)/NUW(I)
      TAWL(I)= TSE(I)*(1.+PR**(1./2.)*(TEM1-1.))
      TAWT(I)= TSE(I)*(1.+PR**(1./3.)*(TEM1-1.))
      TBAR(I)= .5*(TWAL(I)+TSE(I))+.22*PR**(1./3.)*(TTZ-TSE(I))
      MUBAR(I)= MUTZ*SUTHL(I)*TBAR(I)/TTZ
      BB(I) = ME(I)*ATZ/NUTZ*(TSE(I)/TTZ)**(GMP1/(2.*GMM1))
      AA(I)= BB(I)*TSE(I)/TBAR(I)*(MUBAR(I)/MUTZ)**.268
      FF(I)=1.+.1599*ME(I)**2+.60*SW(I)+.2101*SW(I)*ME(I)**2+.0114*ME(I)
     1**4+.0180*SW(I)*ME(I)**4+.1825*SW(I)**2+.0735*SW(I)**2*ME(I)**2
     2+.0073*SW(I)**2*ME(I)**4
  220 CONTINUE
C
C    FINITE DIFFERENCE METHOD USED TO CALCULATE VELOCITY AND MACH NUMBER
C    GRADIENTS ALONG THE SURFACE
C
      DUDS(1) = (UE(2) - UE(1))/(S(2) - S(1))
      DMDS(1) = (ME(2) - ME(1))/(S(2) - S(1))
      DO 230 I=2,NST
      IM = I - 1
      IF (I .EQ. NST) GO TO 230
      IP = I + 1
      DUDS(I) = (UE(IP)-UE(IM))/(S(IP)-S(IM))
      DMDS(I) = (ME(IP)-ME(IM))/(S(IP)-S(IM))
  230 CONTINUE
      DUDS(NST) = (UE(NST) - UE(IM))/(S(NST) - S(IM))
      DMDS(NST) = (ME(NST) - ME(IM))/(S(NST) - S(IM))
  240 DO 250 I=1,NST
  250 DMDL(I)= ARCL*DMDS(I)
C
C  PRINT OTHER CALCULATED PARAMETERS
C
      IF(KPRE.NE.1) GO TO 260
      WRITE(6,1050) (I,AE(I),TSE(I),TWAL(I),TAWL(I),TAWT(I),TBAR(I),
     1I=1,NST)
      WRITE(6,1060) (I,RW(I),SW(I),SUTHL(I),RHSW(I),RHSE(I),HEADW(I),
     1HEADE(I),NUW(I),MUBAR(I),I=1,NST)
  260 IF(KGRAD.NE.1) GO TO 270
      WRITE(6,1070)
      WRITE(6,1080) (I,DUDS(I),DMDS(I),DMDL(I),I=1,NST)
C
C  CHECK FOR IMPROPER INPUT
C
  270 DO 280 I=2,NST
      IF (UE(I).NE.0.) GO TO 280
      ERROR= .TRUE.
      WRITE(6,1090)
      RETURN
  280 CONTINUE
      RETURN
C
C  FORMAT STATEMENTS
C
 1000 FORMAT(1H1///4X,24HPRELIMINARY CALCULATIONS///)
 1010 FORMAT(5X,10HPSZ    =  F12.5/5X,10HTSZ    =  F10.4/5X,10HUZ     =
     1 F11.5//5X,10HASZ    =  F11.4/5X,10HATZ    =  F11.4//5X,10HRHSZ
     2=  G15.7/5X,10HRHTZ   =  G15.7//5X,10HMUSZ   =  G15.7/5X,10HMUTZ
     3 =  G15.7//5X,10HNUSZ   =  G15.7/5X,10HNUTZ   =  G15.7//5X,10HCP
     4   =  F11.4/5X,10HPR     =  F9.5/5X,10HTC     =  G15.7/5X,10HARCL
     5  =  F8.4///)
 1020 FORMAT(/1X,7HSTATION,7X,4HPRES,13X,2HUE,12X,2HME,11X,5HPOPTZ,9X,5H
     1VOVCR/(2X,I3,5X,F12.5,3X,F12.5,4X,F10.6,4X,F10.6,4X,F10.6))
 1030 FORMAT(///1X,7HSTATION,7X,1HX,12X,1HY,12X,1HS,12X,3HXOM,9X,3HYOM,
     19X,3HSOL/(2X,I3,3X,F12.5,1X,F12.5,1X,F12.5,4X,F9.5,3X,F9.5,3X,F9.5
     2))
 1050 FORMAT(///1X,7HSTATION,5X,2HAE,10X,3HTSE,9X,4HTWAL,8X,4HTAWL,8X,4H
     1TAWT,8X,4HTBAR/(2X,I3,4X,F9.3,5(4X,F8.3)))
 1060 FORMAT(///1X,7HSTATION,11X,2HRW,6X,2HSW,4X,5HSUTHL,7X,4HRHSW,12X,4
     1HRHSE,8X,5HHEADW,4X,5HHEADE,9X,3HNUW,12X,5HMUBAR/(2X,I3,3X,F16.1,2
     2X,F4.1,1X,F7.3,2X,G14.6,2X,G14.6,1X,F8.3,1X,F8.3,2X,G14.6,2X,G14.6
     3))
 1070 FORMAT(1H1///21X,17HSURFACE GRADIENTS///)
 1080 FORMAT(1X,7HSTATION,13X,4HDUDS,15X,4HDMDS,15X,4HDMDL/(2X,I3,4X,F18
     1.6,1X,F18.6,1X,F18.6))
 1090 FORMAT(/////10X,83HTHERE IS A STAGNATION POINT AT A STATION OTHER
     1THAN STATION 1.  THIS IS NOT ALLOWED)
      END
$IBFTC LAMNA
      SUBROUTINE LAMNAR

      COMMON/GAMPM/GMP1,GMM1
      COMMON/CTOBL/DUM(2),XMAX,DUM2(5)
      COMMON/C1/GAM,R,PTZ,TTZ,UPMACH,NST,NVP,NTURB,KPVM,KEM,KSMTH,
     1KSPLN,KLE,KATCH,CTHET,DLAM,TLAM,DTURB,TTURB,KPRE,KGRAD,KSDE,KLAM,
     2KMAIN,KPROF,X(100),Y(100),PRES(100),UE(100),ME(100),POPTZ(100),
     3VOVCR(100),TWAL(100),ETA(100)
      COMMON/C2/PSZ,TSZ,UZ,ASZ,ATZ,RHSZ,RHTZ,MUSZ,MUTZ,NUSZ,NUTZ,CP,
     1PR,TC,ARCL
      COMMON/C3/XOM(100),YOM(100),S(100),SOL(100),AE(100),TSE(100),
     1TAWL(100),TAWT(100),TBAR(100),RW(100),SW(100),SUTHL(100),
     2RHSW(100),RHSE(100),HEADW(100),HEADE(100),NUW(100),MUBAR(100),
     3AA(100),BB(100),FF(100),DUDS(100),DMDS(100),DMDL(100)
      COMMON/C4/THET(100),DELSR(100),DELTA(100),FORM(100),
     1FORMI(100),FORMTR(100),RTH(100),RTHI(100),CF(100),
     1TAUW(100),NUSS(100),DTDY(100),HTRAN(100),CRN(100)
      COMMON/C5/SHAPL(100),SHAPK(100),B,NS
      COMMON/C6/FTRAN,FORMS
      COMMON/C7/INST,ITRAN,ISEP
      COMMON/C9/ERROR,TRANS,SEPRN

      DIMENSION CORLN(100),CORML(100),SHEAR(100),DTH(100)
      DIMENSION CCN(20),CRCR(20),CDIF(20),CSHR(20),CCRN(20),CDTH(20)
      DIMENSION STAB(505),CTAB1(505),CTAB2(505)
      REAL MUSZ,NUSZ,MUTZ,NUTZ,ME,NUW,MUBAR,NUSS,NURW,KBAR,INT1,INT2
      LOGICAL ERROR,TRANS,SEPRN
      EXTERNAL FUNCT,INT1,INT2
C
C  READ DATA FOR CORLN(1), RCRIT, DIFF, SHEAR, CRN, AND DTH CURVE FITS
C
      DATA(CCN(I),I=1,6)/-.08178,.06670,-.03143,
     1.00873,.01657,-.01052/
      DATA(CRCR(I),I=1,6)/5.47073,43.6053,227.198,
     1-2067.04,-27172.7,13691.2/
      DATA(CDIF(I),I=1,6)/903.785,26365.0,3.85695E+5,
     11.11044E+6,-4.53853E+7,-7.70276E+7/
      DATA(CSHR(I),I=1,16)/.224488,-1.91539,-9.894,-68.13488,
     1-.001512,-1.4768,-10.52925,-152.2781,-.002406,-.015629,
     1-1.45743,-126.23395,.000752,.005385,.917838,-39.40644/
      DATA(CCRN(I),I=1,16)/2.02056,-19.7211,-24.0495,-1400.002,
     1-.050979,-10.88012,62.4419,-5081.76,-.014343,2.279845,
     1129.7008,-6257.848,.0270567,-1.677051,57.4397,-2552.266/
      DATA(CDTH(I),I=1,16)/8.02829,-4.30978,88.8244,36.4336,
     12.71101,-7.42259,242.293,-16.293, -.16394,-7.61942,286.9795,
     164.11186,-.16758,-3.70289,130.8107,111.3276/
C
C  INITIALIZE PARAMETERS
C
      INST = 0
      ITRAN = 0
      ISEP = 0
      CF(1)= 0.
      TAUW(1)= 0.
      NUSS(1)= 0.
      DTDY(1)= 0.
      HTRAN(1)= 0.
      CRN(1)= 0.
      RTRAN= 0.
      IF (CTHET .GT. 0.) KATCH = 1
C
C  CHECK CONSISTENCY OF INITIAL VALUES
C
      IF (DLAM.GE.0..AND.TLAM.GE.0..AND.DTURB.GE.0..AND.TTURB.GE.0.)
     1GO TO 10
      ERROR = .TRUE.
      WRITE(6,1000)
      RETURN
   10 IF (NTURB.NE.1) GO TO 30
      ITRAN = 1
      IF (DTURB.GT.0..AND.TTURB.GT.0.) GO TO 20
      ERROR = .TRUE.
      WRITE(6,1010)
      RETURN
   20 IF (UE(1).GT.0.) GO TO 240
      ERROR = .TRUE.
      WRITE(6,1020)
      RETURN
C
C  BEGIN CALCULATION IN LAMINAR REGION - CHECK FOR INITIAL VALUES
C  CALCULATE INITIAL CORRELATION NUMBER
C
   30 IF (DLAM.EQ.0..AND.TLAM.EQ.0.) GO TO 70
      IF (UE(1).GT.0.) GO TO 40
      ERROR = .TRUE.
      WRITE(6,1030)
      RETURN
   40 IF (TLAM.EQ.0.) GO TO 50
C  INITIAL MOMENTUM THICKNESS WAS GIVEN
      TEM1 = 1.+GMM1/2.*ME(1)**2
      CORML(1)= -ATZ*TLAM**2/NUTZ/SUTHL(1)/ARCL/TEM1**((3.-GAM)/
     1(2.*GMM1))
      CORLN(1) = CORML(1)*DMDL(1)
      GO TO 90
C  INITIAL DISPLACEMENT THICKNESS WAS GIVEN
   50 IF (ABS(DMDL(1)).GE..0001) GO TO 60
      CORLN(1)= 0.
      TEM1 = 1.+GMM1/2.*ME(1)**2
      FORM(1)= 2.38411*(1.+(2.79-1.78*PR**.5)*((1.+SW(1))*TEM1-1.))+(4.6
     15*PR**(1./3.)-3.65*PR*.5)*PR**.5*(TEM1-1.)
      THET(1)= DLAM/FORM(1)
      CORML(1) = -ATZ*THET(1)**2/NUTZ/SUTHL(1)/ARCL/TEM1**((3.-GAM)/(2.*
     1GMM1))
      GO TO 90
   60 IF(DMDL(1).GT.0.) CALL ROOTB(-1.,0.,DLAM,FUNCT,.5E-5,CORLN(1),SL)
      IF(DMDL(1).LT.0.) CALL ROOTB(0.,.2,DLAM,FUNCT,.5E-5,CORLN(1),SL)
      CORML(1) = CORLN(1)/DMDL(1)
      GO TO 90
C
C  NO INITIAL LAMINAR VALUES GIVEN
C  CALCULATE INITIAL CORRELATION NUMBER
C
C  SHARP LEADING EDGE
   70 IF(KLE.NE.1.AND.ABS(DMDL(1)).GE..0001) GO TO 80
      CORLN(1)= 0.
      CORML(1)= 0.
      GO TO 90
C  STAGNATION POINT
   80 CALL CURVFT(CCN,CORLN(1),SW(1),0,5,0)
      CORML(1)= CORLN(1)/DMDL(1)
      IF (CORML(1).LT.0.) GO TO 90
      ERROR= .TRUE.
      WRITE(6,1040)
      RETURN
C
C  SOLVE LAMINAR DIFFERENTIAL EQUATION
C  CALCULATE CORRELATION NUMBERS ALONG THE SURFACE
C
   90 TEM1 = 1.+GMM1/2.*ME(1)**2
      TEM2 = (3.*GAM-1.)/(2.*GMM1)
      DEL= 0.002*ARCL
      SS= -DEL
      NTAB=1
      CTAB1(1)= CORLN(1)
      CTAB2(1)= CORML(1)
      STAB(1)= 0.
  100 SS= SS+DEL
      SSDEL = SS+DEL
      CALL LGRNGE(S,SW,NST,SS,ANS1)
      CALL LGRNGE(S,ME,NST,SS,ANS2)
      CALL LGRNGE(S,ME,NST,SSDEL,ANS3)
      CALL LGRNGE(S,DMDL,NST,SSDEL,ANS4)
      A1= 0.43631-0.00367*ANS1+0.00481*ANS1**2+0.00651*ANS1**3
      A2= 5.43220+2.25400*ANS1-0.06672*ANS1**2-0.20637*ANS1**3
      A3= 4.51903-10.49775*ANS1-12.71732*ANS1**2-2.95270*ANS1**3
      A4= 19.01831+62.76597*ANS1+115.00986*ANS1**2+62.53113*ANS1**3
      A= A1-A3*CTAB1(NTAB)**2-2.*A4*CTAB1(NTAB)**3
      B= A2+2.*A3*CTAB1(NTAB)+3.*A4*CTAB1(NTAB)**2
C
C    FOR SW = 0.0
      IF ( CTAB1(NTAB).GE.-.1) GO TO 101
      A=.3953
      B=4.739
  101 K1 = 0
      SOL1 = SS/ARCL
      SOL2 = SSDEL/ARCL
      TEM3 = SIMPS1(SOL1,SOL2,INT1,K1)
      IF (TEM3.EQ.0..OR.K1.EQ.0) GO TO 110
      ERROR= .TRUE.
      WRITE(6,1050)
      RETURN
  110 IF (NTAB.GT.1) TEM4= ANS2**(-B)*TEM1**TEM2
      TEM1 = 1.+GMM1/2.*ANS3**2
      TEM5= ANS3**(-B)*TEM1**TEM2
      TEM6= -A*TEM5*TEM3
      IF (NTAB.EQ.1) TEM7=0.
      IF (NTAB.GT.1) TEM7= TEM5/TEM4*CTAB2(NTAB)
      NTAB= NTAB+1
      CTAB2(NTAB)= TEM6+TEM7
      CTAB1(NTAB)= CTAB2(NTAB)*ANS4
      STAB(NTAB)= SSDEL
C
C    WHEN SW IS NOT EQUAL TO 0.0 ,CURVE FIT RANGE ON  CORLN
C    IS FROM -.32 TO .16
C
      IF (CTAB1(NTAB) .GT. .50) GO TO 120
      IF (SS.LT.ARCL) GO TO 100
  120 IF (KSDE.NE.1) GO TO 130
      WRITE(6,1060)
      WRITE(6,1070) (STAB(I),CTAB1(I),I=1,NTAB)
C
C  CALCULATE LAMINAR BOUNDARY LAYER PARAMETERS AT EACH STATION
C
  130 IF (KLAM.NE.1) GO TO 140
      WRITE(6,1080)
  140 I= 0
  150 I= I+1
      IF (I.EQ.NTURB) ITRAN=-1
      IF (S(I).LE.STAB(NTAB)) GO TO 160
      WRITE(6,1090)
  160 IF (KLE .EQ. 1  .AND.  I .EQ. 1) GO TO 151
      CALL LGRNGE(STAB,CTAB1,NTAB,S(I),CORLN(I))
      CALL LGRNGE(STAB,CTAB2,NTAB,S(I),CORML(I))
C  OBTAIN SHEAR, CRN, AND DTH FROM CURVE FITS VS CORLN AND SW
  151 CALL CURVFT(CSHR,SHEAR(I),CORLN(I),SW(I),3,3)
      CALL CURVFT(CCRN,CRN(I),CORLN(I),SW(I),3,3)
      CALL CURVFT(CDTH,DTH(I),CORLN(I),SW(I),3,3)
C
C    FOR SW = 0.0
      IF(CORLN(I).GE. -.1) GO TO 161
      SHEAR(I) = -1.2222*CORLN(I)+.26
      CRN(I) = -58.824 *CORLN(I)-.6765
      DTH(I) = -22.222*CORLN(I)+7.1112
C  CALCULATE OTHER LAMINAR BOUNDARY LAYER PARAMETERS
  161 TEM1 = 1.+GMM1/2.*ME(I)**2
      THET(I)= SQRT(-CORML(I)*NUTZ*SUTHL(I)*ARCL/ATZ*TEM1**((3.-GAM)/
     1(2.*GMM1)))
      FORM(I)= (-1.1138*CORLN(I)+2.38411)*(1.+(2.79-1.78*PR**.5)*((1.+
     1SW(I))*TEM1-1.))+(4.65*PR**(1./3.)-3.65*PR**.5)*PR**.5*(TEM1-1.)
      DELSR(I)= THET(I)*FORM(I)
      RTH(I)= UE(I)*THET(I)/NUW(I)
      FORMI(I)= (FORM(I)-SQRT(PR)*(TEM1-1.))/((1.+SW(I))*TEM1)
      FORMTR(I)= FORMI(I)*(1.+SW(I))
      DELTA(I)= THET(I)*(DTH(I)+(TEM1-1.)*(FORMTR(I)+1.))
      SHAPL(I)= DELTA(I)**2/NUW(I)*DUDS(I)
      IF (I.EQ.1) GO TO 180
      CFRW= 2.*SHEAR(I)*SQRT(-SOL(I)/ME(I)/CORML(I))
      CF(I)= CFRW/SQRT(RW(I))
      TAUW(I)= CF(I)*HEADW(I)
      NURW= CFRW*PR**.3/CRN(I)
      NUSS(I)= NURW*SQRT(RW(I))
      DTDY(I)= NUSS(I)*(TAWL(I)-TWAL(I))/S(I)
      HTRAN(I)= TC*DTDY(I)
      IF (TAUW(I).GT.0.) GO TO 180
      IF (KATCH.NE.0) GO TO 170
      ISEP= I
      SEPRN= .TRUE.
      RETURN
  170 ITRAN= -2
      GO TO 270
  180 IF (I.EQ.1.AND.UE(1).EQ.0.) GO TO 190
      SHAPK(I)= NUTZ*RTH(I)**2*SUTHL(I)**2*(1.+SW(I))**4/ATZ/ME(I)**2/
     1FF(I)/ARCL*DMDL(I)*TEM1**(1./GMM1)
      GO TO 200
  190 SHAPK(I)= 0.07
  200 RTHI(I)= RTH(I)*SUTHL(I)*(1.+SW(I))**2/FF(I)/SQRT(TEM1)
C
C  CALCULATE RCRIT TO CHECK FOR INSTABILITY AND TRANSITION
C
      CALL CURVFT(CRCR,RCRIT,SHAPK(I),0,5,0)
      IF (SHAPK(I) .GT. .07) RCRIT = 8.3163
      RCRIT= EXP(RCRIT)
      IF(INST.NE.0) GO TO 210
C
C  CHECK FOR INSTABILITY
C
      IF(RTHI(I).LT.RCRIT) GO TO 270
      RINS= RTHI(I)
      INST= I
      GO TO 270
C
C  CHECK FOR TRANSITION
C
  210 K1= 0
      NS= I
      TEM= SIMPS1(SOL(INST),SOL(I),INT2,K1)
      IF (TEM.EQ.0..OR.K1.EQ.0) GO TO 220
      ERROR= .TRUE.
      WRITE(6,1100)
      RETURN
  220 KBAR= TEM/(SOL(I)-SOL(INST))
      CALL CURVFT(CDIF,DIFF,KBAR,0,5,0)
      IF (KBAR .GT. .03) DIFF = 44000.*KBAR+700.0
      RTRAN= RINS+DIFF
      IF(RTHI(I).LT.RTRAN) GO TO 270
      IF (I.LT.NTURB) GO TO 270
      ITRAN= -1
      GO TO 270
  230 ITRAN= I
C
C  COMPUTE INITIAL VALUES FOR TURBULENT SOLUTION
C
  240 TRANS= .TRUE.
      IF (DTURB.EQ.0..AND.TTURB.EQ.0.) GO TO 260
      IF (DTURB.GT.0..AND.TTURB.GT.0.) GO TO 250
      ERROR = .TRUE.
      WRITE(6,1110)
      RETURN
  250 THET(ITRAN)= TTURB
      FORM(ITRAN)= DTURB/TTURB
      TEM1 = 1.+GMM1/2.*ME(ITRAN)**2
      FORMI(ITRAN)= (FORM(ITRAN)-PR**(1./3.)*(TEM1-1.))/((1.+SW(ITRAN))
     1*TEM1)
  260 IF (CTHET.GT.0..AND.DTURB.EQ.0..AND.TTURB.EQ.0.) THET(ITRAN)=
     1CTHET*THET(ITRAN)
      THETTR = THET(ITRAN)*(TSE(ITRAN)/TTZ)**(GMP1/(2.*GMM1))
      FTRAN= (ME(ITRAN)*ATZ*THETTR/NUTZ)**1.268
      IF (RTRAN .LE. 0.) GO TO 265
      FORMS= FORMI(ITRAN)-0.59389-0.06591*ALOG(RTRAN)+0.001272*(ALOG(RTR
     1AN))**2
      IF (DTURB.GT.0..AND.TTURB.GT.0.) FORMS=FORMI(ITRAN)
      RETURN
  265 FORMS = 1.4
      RETURN
C
C  PRINT OUTPUT
C
  270 IF (KLAM.NE.1) GO TO 280
      IF(INST.EQ.0 .OR. INST.EQ.I) WRITE(6,1120) I,CORLN(I),SHEAR(I),
     1DTH(I),FORMTR(I),SHAPL(I),RTHI(I),SHAPK(I),RCRIT
      IF(INST.NE.0 .AND. INST.NE.I) WRITE(6,1130) I,CORLN(I),SHEAR(I),
     1DTH(I),FORMTR(I),SHAPL(I),RTHI(I),KBAR,DIFF,RTRAN
      IF (ITRAN.EQ.-2) WRITE(6,1140)
  280 IF(ITRAN.EQ.-1.OR.ITRAN.EQ.-2) GO TO 230
      IF (I.EQ.NST) RETURN
      GO TO 150
C
C  FORMAT STATEMENTS
C
 1000 FORMAT(/////,10X,60HA NEGATIVE INITIAL VALUE HAS BEEN GIVEN. THIS
     1IS NOT ALLOWED)
 1010 FORMAT(/////,10X,75HINITIAL VALUES WERE NOT GIVEN FOR THE TURBULEN
     1T BOUNDARY LAYER AT STATION 1)
 1020 FORMAT(/////,10X,80HINITIAL VALUES WERE GIVEN FOR THE TURBULENT BO
     1UNDARY LAYER AT A STAGNATION POINT)
 1030 FORMAT(/////,10X,94HINITIAL VALUES OTHER THAN ZERO WERE GIVEN FOR
     1THE LAMINAR BOUNDARY LAYER AT A STAGNATION POINT)
 1040 FORMAT(/////,10X,106HFOR THIS INPUT DATA STATION 1 IS ASSUMED TO B
     1E A STAGNATION POINT, SINCE NO INITIAL THICKNESSES ARE GIVEN./
     210X,118HIN THIS CASE PRESSURE SHOULD DECREASE INITIALLY.  EITHER G
     3IVE AN INITIAL VALUE FOR DISPLACEMENT OR MOMENTUM THICKNESS,/
     410X,60HOR BEGIN WITH A SHORT REGION OF FAVORABLE PRESSURE GRADIENT
     5.)
 1050 FORMAT(/////,10X,37HERROR IN COMPUTING INTEGRAL FOR CORLN)
 1060 FORMAT(1H1///7X,50HLAMINAR DIFFERENTIAL EQUATION - SOLUTION FOR CO
     1RLN///5(24H       S       CORLN    )//)
 1070 FORMAT((5(F12.5,2X,F7.4,3X)))
 1080 FORMAT(1H1///1X,59HLAMINAR CALCULATION OF INSTABILITY AND TRANSITI
     1ON LOCATIONS///1X,7HSTATION,2X,5HCORLN,5X,5HSHEAR,5X,3HDTH,6X,6HFO
     2RMTR,4X,5HSHAPL,9X,4HRTHI,6X,5HSHAPK,9X,5HRCRIT,6X,4HKBAR,10X,4HDI
     3FF,9X,5HRTRAN)
 1090 FORMAT(/////,10X,65HLAMINAR SOLUTION HAS PROCEEDED BEYOND THE RANG
     1E WHERE IT IS VALID)
 1100 FORMAT(/////,10X,36HERROR IN COMPUTING INTEGRAL FOR KBAR)
 1110 FORMAT(/////,10X,64HIF INITIAL TURBULENT VALUES ARE GIVEN, THEY BO
     1TH MUST BE NONZERO)
 1120 FORMAT(I4,1X,5F10.4,1X,F12.1,1X,F10.5,1X,F12.1)
 1130 FORMAT(I4,1X,5F10.4,1X,F12.1,24X,F12.5,1X,F12.1,1X,F12.1)
 1140 FORMAT(/////,10X,85HLAMINAR SEPARATION HAS OCCURRED. ASSUMED TO BE
     1 TRANSITION TO TURBULENT BOUNDARY LAYER)
      END
$IBFTC TURBL
      SUBROUTINE TURBLN

      COMMON/GAMPM/GMP1,GMM1
      COMMON/C1/GAM,R,PTZ,TTZ,UPMACH,NST,NVP,NTURB,KPVM,KEM,KSMTH,
     1KSPLN,KLE,KATCH,CTHET,DLAM,TLAM,DTURB,TTURB,KPRE,KGRAD,KSDE,KLAM,
     2KMAIN,KPROF,X(100),Y(100),PRES(100),UE(100),ME(100),POPTZ(100),
     3VOVCR(100),TWAL(100),ETA(100)
      COMMON/C2/PSZ,TSZ,UZ,ASZ,ATZ,RHSZ,RHTZ,MUSZ,MUTZ,NUSZ,NUTZ,CP,
     1PR,TC,ARCL
      COMMON/C3/XOM(100),YOM(100),S(100),SOL(100),AE(100),TSE(100),
     1TAWL(100),TAWT(100),TBAR(100),RW(100),SW(100),SUTHL(100),
     2RHSW(100),RHSE(100),HEADW(100),HEADE(100),NUW(100),MUBAR(100),
     3AA(100),BB(100),FF(100),DUDS(100),DMDS(100),DMDL(100)
      COMMON/C4/THET(100),DELSR(100),DELTA(100),FORM(100),
     1FORMI(100),FORMTR(100),RTH(100),RTHI(100),CF(100),
     1TAUW(100),NUSS(100),DTDY(100),HTRAN(100),CRN(100)
      COMMON/C6/FTRAN,FORMS
      COMMON/C7/INST,ITRAN,ISEP
      COMMON/C8/XTAB(505),YTAB1(505),YTAB2(505),NTAB
      COMMON/C9/ERROR,TRANS,SEPRN

      REAL MUSZ,NUSZ,MUTZ,NUTZ,ME,NUW,MUBAR,NUSS
      LOGICAL ERROR,TRANS,SEPRN
C
C  SOLVE TURBULENT BOUNDARY LAYER DIFFERENTIAL EQUATIONS
C  USING RUNGA-KUTTA
      CALL RUNKUT
      IF (KSDE.NE.1) GO TO 10
      WRITE(6,1000)
      WRITE(6,1010) (XTAB(I),YTAB1(I),YTAB2(I),I=1,NTAB)
   10 DO 5 I=1,NTAB
      IF (YTAB2(I) .LE. 2.8) GO TO 5
      WRITE (6,100)
  100 FORMAT (9H0(TURBLN),5X,63HINCOMPRESSIBLE FORM FACTOR GREATER THAN
     *2.8  -  CASE TERMINATED)
      STOP
    5 CONTINUE
C
C  CALCULATE TURBULENT BOUNDARY LAYER PARAMETERS AT EACH STATION
C
      DO 30 I=ITRAN,NST
      IF (S(I).LE.XTAB(NTAB)) GO TO 20
      ISEP = I-1
      SEPRN= .TRUE.
      RETURN
   20 TEM1 = 1.+GMM1/2.*ME(I)**2
      CALL LGRNGE(XTAB,YTAB1,NTAB,S(I),F)
      THETTR= NUTZ*F**.7886/ME(I)/ATZ
      THET(I) = THETTR*(TTZ/TSE(I))**(GMP1/(2.*GMM1))
      RTH(I)= UE(I)*THET(I)/NUW(I)
      CALL LGRNGE(XTAB,YTAB2,NTAB,S(I),FORMI(I))
      FORMTR(I)= FORMI(I)*(1.+SW(I))
      FORM(I)= FORMTR(I)*TEM1+PR**(1./3.)*(TEM1-1.)
      DELSR(I)= THET(I)*FORM(I)
      POWER= 2.0/(FORMI(I)-1.0)
      IF (FORMI(I).LT.1.02) POWER=100.
      DELTA(I)= (1.+POWER)*DELSR(I)
      CF(I)= 0.246*EXP(-1.561*FORMI(I))*(UE(I)*THET(I)/NUTZ/(TEM1**(1./(
     1GAM-1.))))**(-.268)*TSE(I)/TBAR(I)*(MUBAR(I)/MUTZ)**(.268)
      TAUW(I)= CF(I)*HEADE(I)
      IF (I.EQ.1) GO TO 30
      HTRAN(I)= CF(I)/2./PR**(2./3.)*RHSE(I)*UE(I)*CP*(TAWT(I)-TWAL(I))
      DTDY(I)= HTRAN(I)/TC
      NUSS(I)= S(I)*DTDY(I)/(TAWT(I)-TWAL(I))
      CRN(I)= CF(I)*RW(I)/NUSS(I)
   30 CONTINUE
      RETURN
 1000 FORMAT(1H1///5X,62HTURBULENT DIFFERENTIAL EQUATIONS - SOLUTION FOR
     1  F  AND  FORMI///4(31H      S         F       FORMI  )//)
 1010 FORMAT((4(F10.5,2X,F8.1,2X,F7.4,2X)))
      END
$IBFTC PROFI
      SUBROUTINE PROFIL

      COMMON/CODE/KODE,DELSRL,DELSRU
      COMMON/CTOBL/DUM(2),XMAX,DUM2(5)
      COMMON/LNK2/DUM3(3),CSTAR,DUM35(2),GSTARI,DUM4(3)
      COMMON/C1/GAM,R,PTZ,TTZ,UPMACH,NST,NVP,NTURB,KPVM,KEM,KSMTH,
     1KSPLN,KLE,KATCH,CTHET,DLAM,TLAM,DTURB,TTURB,KPRE,KGRAD,KSDE,KLAM,
     2KMAIN,KPROF,X(100),Y(100),PRES(100),UE(100),ME(100),POPTZ(100),
     3VOVCR(100),TWAL(100),ETA(100)
      COMMON/C3/XOM(100),YOM(100),S(100),SOL(100),AE(100),TSE(100),
     1TAWL(100),TAWT(100),TBAR(100),RW(100),SW(100),SUTHL(100),
     2RHSW(100),RHSE(100),HEADW(100),HEADE(100),NUW(100),MUBAR(100),
     3AA(100),BB(100),FF(100),DUDS(100),DMDS(100),DMDL(100)
      COMMON/C4/THET(100),DELSR(100),DELTA(100),FORM(100),
     1FORMI(100),FORMTR(100),RTH(100),RTHI(100),CF(100),
     1TAUW(100),NUSS(100),DTDY(100),HTRAN(100),CRN(100)
      COMMON/C5/SHAPL(100),SHAPK(100),B,NS
      COMMON/C7/INST,ITRAN,ISEP
      REAL ME,NUSS
C
C  PRINT LOCATIONS OF INSTABILITY, TRANSITION, AND SEPARATION
C
      IF (KMAIN.NE.1) GO TO 60
      WRITE(6,1000)
      IF(INST.EQ.0) GO TO 10
      WRITE(6,1010) INST
      GO TO 20
   10 WRITE(6,1020)
   20 IF (ITRAN.LE.1) GO TO 30
      WRITE(6,1030) ITRAN
      GO TO 40
   30 WRITE(6,1040)
   40 IF(ISEP.EQ.0) GO TO 50
      WRITE(6,1050) ISEP
      GO TO 60
   50 WRITE(6,1060)
C
C  PRINT LOCATIONS OF LAMINAR AND TURBULENT BOUNDARY LAYERS
C
   60 IEND = ITRAN-1
      IF (IEND.EQ.-1.OR.IEND.EQ.0) IEND=ISEP
      IF (IEND.EQ.0) IEND=NST
      IF (KMAIN.NE.1) GO TO 70
      IF (ITRAN.EQ.1) WRITE(6,1070)
      IF (ITRAN.NE.1) WRITE(6,1080) IEND
      IF (ITRAN.EQ.0) WRITE(6,1090)
      IF (ITRAN.EQ.1) WRITE(6,1100) ITRAN,IEND
   70 IF (ITRAN.LE.1) GO TO 80
      IEND = ISEP
      IF(IEND.EQ.0) IEND=NST
      IF (KMAIN.NE.1) GO TO 90
      WRITE(6,1100) ITRAN,IEND

   80 IF (KODE .EQ. 0) DELSRU = DELSR(IEND)/XMAX
      IF (KODE .EQ. 1) DELSRL = DELSR(IEND)/XMAX

      WRITE (6,2000)
 2000 FORMAT (1H1,35X,29HDIMENSIONED BLADE COORDINATES,36X,24HANGLES OF
     *ROTATION (DEG)//37X,11HUNCORRECTED,14X,9HCORRECTED,15X,10HTRANSLAT
     $ED//1X,7HSTATION,12X,1HX,23X,1HY,21X,5HYCORR,18X,5HYTRAN)
      SIGN = +1.
      IF (KODE .EQ. 0) SIGN = -1.
      GSTAR = GSTARI*XMAX/CSTAR
      DO 170 I=1,IEND
      YCORR = Y(I) + SIGN*ABS(DELSR(I)/COS(ETA(I)))
      YTRAN = YCORR - GSTAR
      ETAD = 57.2957796*ETA(I)
      WRITE (6,2001) I,X(I),Y(I),YCORR,YTRAN,ETAD
 2001 FORMAT (2X,I3,9X,4(F9.5,15X),F9.4)
  170 CONTINUE
C
C  PRINT CALCULATED BOUNDARY LAYER PARAMETERS
C
      IF (KMAIN .NE. 1) GO TO 90
      WRITE(6,1110)
      WRITE(6,1120) (I,X(I),S(I),DELSR(I),THET(I),DELTA(I),FORM(I),
     1FORMI(I),I=1,IEND)
      WRITE(6,1130)
      WRITE(6,1140) (I,CF(I),TAUW(I),RTH(I),DTDY(I),NUSS(I),HTRAN(I),
     1CRN(I),I=1,IEND)
C
C  COMPUTE BOUNDS ON VELOCITY PROFILES
C
   90 IF (KPROF.NE.1) RETURN
      WRITE(6,1150)
      IF(ITRAN.NE.0) GO TO 100
      IL1= 2
      IL2= IEND
      IT1= 0
      IT2= 0
      GO TO 110
  100 IL1= 2
      IL2= ITRAN-1
      IT1= ITRAN
      IT2= IEND
      IF (IT1.EQ.1) IT1=2
C
C  CALCULATE AND PRINT LAMINAR BOUNDARY LAYER VELOCITY PROFILES
C
  110 NVP1= NVP+1
      IF (IL2.LT.IL1) GO TO 140
      DO 130 I=IL1,IL2
      WRITE(6,1160) I
      AAA= 2.+SHAPL(I)/6.
      BBB= -.5*SHAPL(I)
      CCC= -2.+.5*SHAPL(I)
      DDD= 1.-SHAPL(I)/6.
      DEL= DELTA(I)/FLOAT(NVP)
      YP= -DEL
      DO 120 J=1,NVP1
      YP= YP+DEL
      ETAA = YP/DELTA(I)
      YXMAX= YP/X(NST)
      UUE = (((DDD*ETAA+CCC)*ETAA+BBB)*ETAA+AAA)*ETAA
      U= UUE*UE(I)
  120 WRITE (6,1180) ETAA,YP,YXMAX,U,UUE
  130 CONTINUE
C
C  CALCULATE AND PRINT TURBULENT BOUNDARY LAYER VELOCITY PROFILES
C
  140 IF(IT1.EQ.0) RETURN
      DO 160 I=IT1,IT2
      POWER= DELTA(I)/DELSR(I)-1.
      WRITE(6,1170) I,POWER
      DEL= DELTA(I)/FLOAT(NVP)
      YP= -DEL
      DO 150 J=1,NVP1
      YP= YP+DEL
      ETAA = YP/DELTA(I)
      YXMAX = YP/XMAX
      UUE = ETAA**(1./POWER)
      U= UUE*UE(I)
  150 WRITE (6,1180) ETAA,YP,YXMAX,U,UUE
  160 CONTINUE
      RETURN
C
C  FORMAT STATEMENTS
C
 1000 FORMAT(1H1///1X,36HPRINCIPAL BOUNDARY LAYER INFORMATION///)
 1010 FORMAT (/10X,31HINSTABILITY OCCURS AT STATION  ,I3)
 1020 FORMAT (/10X,26HINSTABILITY DOES NOT OCCUR)
 1030 FORMAT (/10X,30HTRANSITION OCCURS AT STATION  ,I3)
 1040 FORMAT (/10X,25HTRANSITION DOES NOT OCCUR)
 1050 FORMAT (/10X,30HSEPARATION OCCURS AT STATION  ,I3)
 1060 FORMAT (/10X,25HSEPARATION DOES NOT OCCUR)
 1070 FORMAT (/10X,37HLAMINAR BOUNDARY LAYER DOES NOT OCCUR)
 1080 FORMAT (/10X,42HLAMINAR BOUNDARY LAYER - STATIONS  1  TO  ,I3)
 1090 FORMAT (/10X,39HTURBULENT BOUNDARY LAYER DOES NOT OCCUR///)
 1100 FORMAT (/10X,35HTURBULENT BOUNDARY LAYER - STATIONS,2X,
     1I3,6H  TO  ,I3///)
 1110 FORMAT(/1X,7HSTATION,8X,1HX,12X,1HS,12X,5HDELSR,10X,4HTHET,11X,
     15HDELTA,11X,4HFORM,10X,5HFORMI)
 1120 FORMAT(2X,I3,3X,2F13.6,F14.6,1X,F14.6,1X,F14.6,1X,2F14.4)
 1130 FORMAT(///1X,7HSTATION,6X,2HCF,13X,4HTAUW,11X,3HRTH,14X,4HDTDY,
     113X,4HNUSS,10X,5HHTRAN,12X,3HCRN)
 1140 FORMAT(I5,F14.5,2X,F14.5,1X,F12.1,5X,F14.2,2X,F14.2,1X,
     1F14.4,2X,F13.3)
 1150 FORMAT(1H1///1X,17HVELOCITY PROFILES///)
 1160 FORMAT (/1X,7HSTATION,1X,I5,2X,7HPROFILE/3X,7HY/DELTA,9X,
     11HY,12X,6HY/XMAX,10X,1HU,12X,4HU/UE)
 1170 FORMAT (/1X,7HSTATION,1X,I5,2X,7HPROFILE,28X,2HN=,1X,F6.2/3X,7HY/D
     1ELTA,9X,1HY,12X,6HY/XMAX,10X,1HU,12X,4HU/UE)
 1180 FORMAT(1X,F8.4,2X,2G15.6,2X,F9.2,6X,F8.4)
      END
$IBFTC RUNKT
      SUBROUTINE RUNKUT
C
C  RUNKUT SOLVES SIMULTANEOUS FIRST ORDER INITIAL VALUE
C  ORDINARY DIFFERENTIAL EQUATIONS
C
      COMMON/C1/GAM,R,PTZ,TTZ,UPMACH,NST,NVP,NTURB,KPVM,KEM,KSMTH,
     1KSPLN,KLE,KATCH,CTHET,DLAM,TLAM,DTURB,TTURB,KPRE,KGRAD,KSDE,KLAM,
     2KMAIN,KPROF,X(100),Y(100),PRES(100),UE(100),ME(100),POPTZ(100),
     3VOVCR(100),TWAL(100),ETA(100)
      COMMON/C3/XOM(100),YOM(100),S(100),SOL(100),AE(100),TSE(100),
     1TAWL(100),TAWT(100),TBAR(100),RW(100),SW(100),SUTHL(100),
     2RHSW(100),RHSE(100),HEADW(100),HEADE(100),NUW(100),MUBAR(100),
     3AA(100),BB(100),FF(100),DUDS(100),DMDS(100),DMDL(100)
      COMMON/C6/FTRAN,FORMS
      COMMON/C7/INST,ITRAN,ISEP
      COMMON/C8/XTAB(505),YTAB1(505),YTAB2(505),NTAB
      DIMENSION YY(2),RY(2),YINC(2),DOT(2),RUK(2,4)
      DOUBLE PRECISION XX,RX,YY,RY,RUK,DEL,DOT,
     1TEM1,TEM2,TEM3,TEM4,TEM5,TEM6
      REAL ME,NUW,MUBAR
C
C  SET DEL SPACING AND STORE INITIAL VALUES
C
      DEL= 0.002*S(NST)
   10 YY(1)=FTRAN
      YY(2)= FORMS
      XX= S(ITRAN)
      NV=2
      NTAB = 1
      YTAB1(1)= YY(1)
      YTAB2(1)= YY(2)
      XTAB(1)= XX
C
C  SOLVE FOR YY(1) AND YY(2) AT NEXT XX INCREMENT
C
C  SAVE PREVIOUS YY(1) AND YY(2)
   20 DO 30 J=1,NV
   30 RY(J)= YY(J)
      RX= XX
C
C  CALCULATE NEW YY(1) AND YY(2)
C
      DO 90 L=1,4
C  PUT DIFFERENTIAL EQUATIONS IN THE FORM OF
C  FIRST DERIVATIVE = REMAINDER OF EQUATION
      CALL LGRNGE(S,ME,NST,XX,ANS1)
      CALL LGRNGE(S,SW,NST,XX,ANS2)
      CALL LGRNGE(S,AA,NST,XX,ANS3)
      CALL LGRNGE(S,BB,NST,XX,ANS4)
      CALL LGRNGE(S,DMDS,NST,XX,ANS5)
      CALL LGRNGE(S,TBAR,NST,XX,ANS6)
      TEM1= 1.+(1.+ANS2)*YY(2)
      TEM2= .123*EXP(-1.561*YY(2))*ANS3
      DOT(1)= 1.268*(-YY(1)/ANS1*ANS5*TEM1+TEM2)
      TEM3= YY(2)*(YY(2)+1.)**2*(YY(2)-1.)
      TEM4= 1.+ANS2*(YY(2)*YY(2)+4.*YY(2)-1.)/((YY(2)+1.)*(YY(2)+3.))
      TEM5= (YY(2)*YY(2)-1.)*YY(2)/YY(1)*(.123*EXP(-1.561*YY(2))*ANS3)
      TEM6= (YY(2)*YY(2)-1.)/YY(1)**(.7886)*(.011*(YY(2)+1.)*(YY(2)-1.)
     1**2/YY(2)**2*TTZ/ANS6)*ANS4
      DOT(2)= -ANS5*.5/ANS1*TEM3*TEM4+TEM5-TEM6
C  APPLY THE RUNGA-KUTTA SCHEME
      DO 40 J=1,NV
   40 RUK(J,L) = DEL*DOT(J)
      GO TO (50,50,70,90), L
   50 DO 60 J=1,NV
   60 YY(J)= RY(J)+RUK(J,L)/2.
      XX= RX+DEL/2.
      GO TO 90
   70 DO 80 J=1,NV
   80 YY(J)= RY(J)+RUK(J,L)
      XX= RX+DEL
   90 CONTINUE
C  INCREMENT THE DEPENDENT VARIABLES TO OBTAIN NEW YY(1) AND YY(2)
      DO 100 J=1,NV
      YINC(J) = (RUK(J,1)+2.*RUK(J,2)+2.*RUK(J,3)+RUK(J,4))/6.
  100 YY(J)= RY(J)+YINC(J)
C
C  STORE NEW COMPUTED VALUES IN A TABLE
C
      NTAB = NTAB+1
      YTAB1(NTAB)= YY(1)
      YTAB2(NTAB)= YY(2)
      XTAB(NTAB)= XX
      IF (YTAB2(NTAB) .GT. 2.8) RETURN
      IF (XX.LT.S(NST)) GO TO 20
      RETURN
      END
$IBFTC SPLIN
      SUBROUTINE SPLINE(X,Y,N,DYDX,D2YDX2)
C
C  SPLINE FITS A SPLINE CURVE TO  X  AND  Y
C  AND CALCULATES FIRST AND SECOND DERIVATIVES AT THE SPLINE POINTS
C  END POINT SECOND DERIVATIVES EQUAL THOSE AT ADJACENT POINTS
C
      DIMENSION X(N),Y(N),DYDX(N),D2YDX2(N)
      DIMENSION G(100),H(100)
      G(1)= -1.
      H(1)= 0.
      N1= N-1
      IF (N1.LT.2) GO TO 20
      DO 10 I=2,N1
      A= (X(I)-X(I-1))/6.
      B= (X(I+1)-X(I))/6.
      C= 2.*(A+B)-A*G(I-1)
      D= (Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1))/(X(I)-X(I-1))
      G(I)= B/C
   10 H(I)= (D-A*H(I-1))/C
   20 D2YDX2(N)= H(N1)/(1.+G(N1))
      DO 30 I=2,N
      K= N+1-I
   30 D2YDX2(K)= H(K)-G(K)*D2YDX2(K+1)
      DYDX(1)= (X(1)-X(2))/6.*(2.*D2YDX2(1)+D2YDX2(2))+(Y(2)-Y(1))/(X(2)
     1-X(1))
      DO 40 I=2,N
   40 DYDX(I)= (X(I)-X(I-1))/6.*(2.*D2YDX2(I)+D2YDX2(I-1))+(Y(I)-Y(I-1))
     1/(X(I)-X(I-1))
      RETURN
      END
$IBFTC LGRNGV
      SUBROUTINE LGRNGE (X,Y,N,ARG,ANS)

C    LINEAR INTERPOLATION (REPLACES 4-POINT LAGRANGE)

      DIMENSION X(N),Y(N)

      IF (ARG .GE. 0.80*X(1)) GO TO 2
      WRITE (6,100) ARG,X(1),X(N)
  100 FORMAT (9H0(LGRNGE),5X,10HABSCISSA =,E13.5,27H IS OUT OF RANGE  -
     * X(1) =,E13.5,5X,6HX(N) =,E13.5/15X,27HEXTRAPOLATED VALUE RETURNED
     *)
    2 NM = N - 1
      DO 3 I=2,NM
      IF (ARG .GT. X(I)) GO TO 3
      M = I
      GO TO 4
    3 CONTINUE
      IF (ARG .LE. 1.20*X(N)) GO TO 5
      WRITE (6,100) ARG,X(1),X(N)
    5 M = N
    4 ANS = Y(M) + (Y(M) - Y(M-1))/(X(M) - X(M-1))*(ARG - X(M))
      RETURN

      END
$IBFTC SIMP
      FUNCTION SIMPS1(X1,X2,FUNC,KSIG)
      DIMENSION V(200),H(200),A(200),B(200),C(200),P(200),E(200)
      LOGICAL SPILL
      DOUBLE PRECISION ANS,Q
      DATA TWO,THREE,FOUR,THIRTY/2.0,3.0,4.0,30.0/
      DATA T,NMAX,NSIG/3.0E-5,200,1/
C  INITIALIZE FIRST ELEMENTS OF ARRAYS.
      V=X1
      H=(X2-V)/TWO
      A=FUNC(V)
      B=FUNC(V+H)
      C=FUNC(X2)
      P=H*(A+FOUR*B+C)
      E=P
      ANS=P
      N=1
      FRAC=T
      SPILL=.FALSE.
   10 TEST=ABS(FRAC*ANS)
      K=N
      DO 30 I=1,K
C  TEST MAGNITUDE OF 4TH ORDER ERROR IN THIS INTERVAL.
      IF (ABS(E(I)).LE.TEST) GO TO 30
      IF (N.LT.NMAX) GO TO 20
C  GO TO FINISH IF STORAGE IS FILLED UP.
      SPILL=.TRUE.
      KSIG=KSIG+NSIG
      GO TO 40
C  SUBDIVIDE INTERVAL AGAIN TO REDUCE 4TH ORDER ERROR.
   20 N=N+1
      V(N)=V(I)+H(I)
      H(N)=H(I)/TWO
      A(N)=B(I)
      B(N)=FUNC(V(N)+H(N))
      C(N)=C(I)
      P(N)=H(N)*(A(N)+FOUR*B(N)+C(N))
      H(I)=H(N)
      B(I)=FUNC(V(I)+H(I))
      C(I)=A(N)
      Q=P(I)
      P(I)=H(I)*(A(I)+FOUR*B(I)+C(I))
      Q=P(I)+P(N)-Q
      ANS=ANS+Q
      E(I)=Q
      E(N)=Q
   30 CONTINUE
C  TEST ALL INTERVALS AGAIN IF ANY WERE SUBDIVIDED THE LAST TIME.
      IF (N.GT.K) GO TO 10
   40 Q=0.0
      DO 50 I=1,N
   50 Q=Q+E(I)
C  TIGHTEN ERROR LIMIT IF TOTAL ACCUMULATED ERROR TOO LARGE.
      IF (ABS(Q/T).LE.ABS(ANS).OR.SPILL) GO TO 60
      FRAC=FRAC/TWO
      GO TO 10
C  FINISH CALCULATION.
   60 SIMPS1=(ANS+Q/THIRTY)/THREE
      RETURN
      END
$IBFTC CURVF
      SUBROUTINE CURVFT(COEF,ANS,X,Y,NX,NY)
C
C  EVALUATE THE POLYNOMIAL FUNCTION, ANS=F(X,Y), USING COEFFICIENTS, COEF
C
      DIMENSION COEF(20)
      NX1 = NX+1
      NY1 = NY+1
      ANS = COEF(1)
      IF (X.EQ..0.AND.Y.EQ..0) RETURN
      IF (Y.EQ..0) GO TO 10
      IF (X.EQ..0) GO TO 30
      GO TO 50
   10 DO 20 I=2,NX1
   20 ANS = ANS+COEF(I)*X**(I-1)
      RETURN
   30 DO 40 I=2,NY1
      K = (I-1)*NX1+1
   40 ANS = ANS+COEF(K)*Y**(I-1)
      RETURN
   50 ANS = .0
      DO 60 I=1,NY1
      DO 60 J=1,NX1
      K = (I-1)*NX1+J
   60 ANS = ANS+COEF(K)*Y**(I-1)*X**(J-1)
      RETURN
      END
$IBFTC AFTMIX
      SUBROUTINE AFMIX (ALPH1,DELS,DELP,THETS,THETP,TE,SP,XMFS1)

      DIMENSION XXX(100),YYY(100),SSS(100)

      COMMON/GAMPM/GMP1,GMM1
      COMMON/C1/GAM,R,PTZ,TTZ,UPMACH,NST,NVP,NTURB,KPVM,KEM,KSMTH,
     1KSPLN,KLE,KATCH,CTHET,DLAM,TLAM,DTURB,TTURB,KPRE,KGRAD,KSDE,KLAM,
     2KMAIN,KPROF,X(100),Y(100),PRES(100),UE(100),ME(100),POPTZ(100),
     3VOVCR(100),TWAL(100),ETA(100)
      COMMON/C3/XOM(100),YOM(100),S(100),SOL(100),AE(100),TSE(100),
     1TAWL(100),TAWT(100),TBAR(100),RW(100),SW(100),SUTHL(100),
     2RHSW(100),RHSE(100),HEADW(100),HEADE(100),NUW(100),MUBAR(100),
     3AA(100),BB(100),FF(100),DUDS(100),DMDS(100),DMDL(100)
      COMMON/C4/THET(100),DELSR(100),DELTA(100),FORM(100),
     1FORMI(100),FORMTR(100),RTH(100),RTHI(100),CF(100),
     1TAUW(100),NUSS(100),DTDY(100),HTRAN(100),CRN(100)
      EQUIVALENCE (X,XXX),(Y,YYY),(S,SSS),(NST,N)
      REAL ME

      KODE = 0
      WRITE(6,5)
5     FORMAT(1H1,20X,22HAFTERMIXING PROPERTIES)
      XALPH1 = ALPH1
      ALPH1 = ALPH1*0.017453
    2 VVCR1 = SQRT((GMP1/2.*XMFS1**2)/(1.+GMM1/2.*XMFS1**2))
      XX=SP*COS(ALPH1)
      DELSRT= (DELS+DELP)/XX
      THETA= (THETS+THETP)/XX
      DTE = TE/XX
      AFS1 = GMM1/GMP1*VVCR1**2
      A = 1.0-DELSRT-DTE-THETA
      A1= 1.0-DELSRT-DTE
      IF (A  .LE. 0.0) GO TO 16
      IF (KODE .EQ. 0) WRITE (6,1)
    1 FORMAT (1H0,10X,39HROTOR WITH NO BOUNDARY LAYER CORRECTION)
      C = ((1.-AFS1)*GMP1/(2.*GAM)+(COS(ALPH1))**2*
     1A*VVCR1**2)/(COS(ALPH1)*A1*VVCR1)
      D = VVCR1*SIN(ALPH1)*A/A1
      VXVCR2 = GAM*C/GMP1-SQRT((GAM*C/GMP1)**2
     1-1.+GMM1/GMP1*D**2)
      GO TO 4
3     KODE=2
      WRITE(6,10)
10    FORMAT(1H0,10X,19HSUPERSONIC SOLUTION)
      VXVCR2 = GAM*C/GMP1+SQRT((GAM*C/GMP1)**2
     1-1.+GMM1/GMP1*D**2)
    4 DENCH = GMM1/GMP1*(D**2+VXVCR2**2)
      IF (DENCH.GE. 1.0) RETURN
      DENSR2 = (1.-GMM1/GMP1*(D**2+VXVCR2**2))
     1**(1./GMM1)
      DENSR1 = 1./((1.+GMM1/2.*XMFS1**2)**(1./GMM1))
      PR2 = DENSR2**GAM
      PT2PT0=(DENSR1*VVCR1*COS(ALPH1)*A1)/(DENSR2*VXVCR2)
      PT0P2 = 1.0/(PT2PT0*PR2)
      EBAR2 = ((PT2PT0)**(-GMM1/GAM)-1.)/(PT0P2**
     1(GMM1/GAM)-1.)
      ETAN = 1.0-EBAR2
      VVCR2 = SQRT(D**2+VXVCR2**2)
      XM2 = SQRT(((2./GMP1)*VVCR2**2)/(1.-(GMM1/GMP1)
     1*VVCR2**2))
      T2TT0 = 1.-GMM1/GMP1*VVCR2**2
      ALPH2 = ATAN(D/VXVCR2)
      XMX1 = XMFS1*COS(ALPH1)
      XMX2 = XM2*COS(ALPH2)
      ALPH2 = ALPH2*57.2958
      WRITE(6,6) XMFS1,SP,TE,XM2,VVCR1,XMX1,XMX2
6     FORMAT(1H0,7HXMFS1 =,F6.4,2X,9HSPACING =,F7.6,2X,4HTE =,F7.5,2X,
     15HXM2 =,F6.4,2X,8HV/VCR1 =,F6.3,2X,6HXMX1 =,F6.3,2X,6HXMX2 =,F6.3)
      WRITE (6,7) XALPH1,ALPH2,PT2PT0,PT0P2,T2TT0,VVCR2,EBAR2,ETAN
7     FORMAT(1H0,6HALPH1=,F7.3,2X,6HALPH2=,F7.3,2X,8HPT2/PT0=,F7.4,2X,
     17HPT0/P2=,F9.3,2X,7HT2/TT0=,F7.4,2X,7HV/VCR2=,F6.3,2X,6HEBAR2=,
     2F7.5,2X,6HETA-N=,F6.4)
      IF (KODE .EQ. 2 ) RETURN
      IF( KODE.EQ.1 ) GO TO 3
   16 SP = SP + (DELS + DELP)/COS(ALPH1)
      WRITE(6,20)
   20 FORMAT (1H0,10X,36HROTOR WITH BOUNDARY LAYER CORRECTION)
      KODE = 1
      GO TO 2
      END
$IBFTC FUNC
      SUBROUTINE FUNCT(XX,FX,DFX,INF)

      COMMON/GAMPM/GMP1,GMM1
      COMMON/C1/GAM,R,PTZ,TTZ,UPMACH,NST,NVP,NTURB,KPVM,KEM,KSMTH,
     1KSPLN,KLE,KATCH,CTHET,DLAM,TLAM,DTURB,TTURB,KPRE,KGRAD,KSDE,KLAM,
     2KMAIN,KPROF,X(100),Y(100),PRES(100),UE(100),ME(100),POPTZ(100),
     3VOVCR(100),TWAL(100),ETA(100)
      COMMON/C2/PSZ,TSZ,UZ,ASZ,ATZ,RHSZ,RHTZ,MUSZ,MUTZ,NUSZ,NUTZ,CP,
     1PR,TC,ARCL
      COMMON/C3/XOM(100),YOM(100),S(100),SOL(100),AE(100),TSE(100),
     1TAWL(100),TAWT(100),TBAR(100),RW(100),SW(100),SUTHL(100),
     2RHSW(100),RHSE(100),HEADW(100),HEADE(100),NUW(100),MUBAR(100),
     3AA(100),BB(100),FF(100),DUDS(100),DMDS(100),DMDL(100)

      REAL MUSZ,NUSZ,MUTZ,NUTZ,ME,NUW,MUBAR

      INF = 0
      B1 = 1.+GMM1/2.*ME(1)**2
      B2= 1.+(2.79-1.78*PR**.5)*((1.+SW(1))*B1-1.)
      B3 = -NUTZ*SUTHL(1)*ARCL/ATZ/DMDL(1)*B1**((3.-GAM)/(2.*GMM1))
      B4= -1.1138*B2
      B5= 2.38411*B2+(4.65*PR**(1./3.)-3.65*PR**.5)*PR**.5*(B1-1.)
      FX= (B3*XX)**.5*(B4*XX+B5)
      IF (XX.EQ.0.) GO TO 10
      DFX= .5*(B3*XX)**(-.5)*B3*(B4*XX+B5)+B4*(B3*XX)**.5
      RETURN
   10 INF = 1
      DFX = 1.E10
      RETURN
      END
$IBFTC ROTB
      SUBROUTINE ROOTB (A,B,Y,FUNCT,TOLERY,X,DFX)
C
C  ROOT FINDS A ROOT FOR (FUNCT-Y) IN THE INTERVAL (A,B)
C
      X1= A
      X2= B
      CALL FUNCT(X1,FX1,DFX,INF)
   10 DO 30 I=1,20
      X= (X1+X2)/2.
      CALL FUNCT(X,FX,DFX,INF)
      IF ((FX1-Y)*(FX-Y).GT.0.) GO TO 20
      X2= X
      GO TO 30
   20 X1= X
      FX1= FX
   30 CONTINUE
      IF (ABS(Y-FX).LT.TOLERY) RETURN
      WRITE(6,1000) A,B,Y
      STOP
 1000 FORMAT(////4X,49HROOT HAS FAILED TO CONVERGE IN THE GIVEN INTERVAL
     1/4X,3HA =,G14.6,10X,3HB =,G14.6,10X,3HY =,G14.6)
      END
$IBFTC INTG1
      REAL FUNCTION INT1(XX)

      COMMON/GAMPM/GMP1,GMM1
      COMMON/C1/GAM,R,PTZ,TTZ,UPMACH,NST,NVP,NTURB,KPVM,KEM,KSMTH,
     1KSPLN,KLE,KATCH,CTHET,DLAM,TLAM,DTURB,TTURB,KPRE,KGRAD,KSDE,KLAM,
     2KMAIN,KPROF,X(100),Y(100),PRES(100),UE(100),ME(100),POPTZ(100),
     3VOVCR(100),TWAL(100),ETA(100)
      COMMON/C3/XOM(100),YOM(100),S(100),SOL(100),AE(100),TSE(100),
     1TAWL(100),TAWT(100),TBAR(100),RW(100),SW(100),SUTHL(100),
     2RHSW(100),RHSE(100),HEADW(100),HEADE(100),NUW(100),MUBAR(100),
     3AA(100),BB(100),FF(100),DUDS(100),DMDS(100),DMDL(100)
      COMMON/C5/SHAPL(100),SHAPK(100),B,NS

      REAL ME,NUW,MUBAR,INT1

      CALL LGRNGE(SOL,ME,NST,XX,ANS)
      INT1 = ANS**(B-1.)/((1.+GMM1/2.*ANS**2)**
     1((3.*GAM-1.)/(2.*GMM1)))
      RETURN
      END
$IBFTC INTG2
      REAL FUNCTION INT2(XX)

      COMMON/C1/GAM,R,PTZ,TTZ,UPMACH,NST,NVP,NTURB,KPVM,KEM,KSMTH,
     1KSPLN,KLE,KATCH,CTHET,DLAM,TLAM,DTURB,TTURB,KPRE,KGRAD,KSDE,KLAM,
     2KMAIN,KPROF,X(100),Y(100),PRES(100),UE(100),ME(100),POPTZ(100),
     3VOVCR(100),TWAL(100),ETA(100)
      COMMON/C3/XOM(100),YOM(100),S(100),SOL(100),AE(100),TSE(100),
     1TAWL(100),TAWT(100),TBAR(100),RW(100),SW(100),SUTHL(100),
     2RHSW(100),RHSE(100),HEADW(100),HEADE(100),NUW(100),MUBAR(100),
     3AA(100),BB(100),FF(100),DUDS(100),DMDS(100),DMDL(100)
      COMMON/C5/SHAPL(100),SHAPK(100),B,NS
      REAL ME,NUW,MUBAR,INT2

      IF (NS.LT.4) GO TO 10
      CALL LGRNGE(SOL,SHAPK,NS,XX,INT2)
      RETURN
   10 DO 20 J=2,NS
      IF (SOL(J).LT.XX) GO TO 20
      INT2= SHAPK(J-1)+(SHAPK(J)-SHAPK(J-1))*(XX-SOL(J-1))/(SOL(J)-SOL(J
     1-1))
      RETURN
   20 CONTINUE
      RETURN
      END
